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The  model  developed  solves  the  multiperiod  multireservoir  water 
resources  problem  with  stochastic  inflows.  Of  unique  importance  is 
the  development  of  a  generalized  network  model  which  solves 
nonlinear  nonseparable  quadratic  problems.  Quadratic  functions  are 
used  to  measure  the  future  value  of  water  to  the  system.  The 
nonseparable  form  stems  from  the  realization  that  interaction 
exists  between  the  benefits  to  be  gained  from  a  multireservoir 
system.  Historically  this  interactive  nature  has  been  ignored  due 
to  the  computational  difficulty  of  measuring  and  solving  such 
relationships.  Also  developed  is  a  3tocha3tic  dynamic  programming 
approach  which  utilizes  the  results  of  the  network  optimization  as 
data  for  a  least  squares  regression  analysis.  A  quadratic  function 
is  fit  to  this  data  and  is  used  to  represent  the  future  value  of 
water  to  the  system  for  the  next  period  in  the  dynamic  programming 
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approach.  This  functional  representation  of  the  future  value  of 
water  replaces  the  standard  discrete  matrix  representation  of 
dynamic  programming  and  greatly  reduces  the  dimensionality  problems 
associated  with  the  dynamic  programming  approach.  In  the  end,  this 
work  represents  a  rare  combination  of  generalized-nonlinear  network 
flow  programming,  stochastic  dynamic  programming  and  regression 
analysis. ^Example  problems  are  included  along  with  an  application 
to  a  four  reservoir  model  of  the  Guadalupe  River  Basin  in  Texas. 
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CHAPTER  I 
1  .  Introduction 


i . 1  General 

During  the  1970’ a  significant  advances  were  made  on 
computational  techniques  for  determining  optimum  solutions  for 
network  flow  problems.  It  is  now  possible  to  solve  problems  of 
tens  of  thousands  of  variables  using  only  seconds  of  time  on  large 
modern  computers.  These  advances,  along  with  their  historical 
precursors  are  described  in  several  recent  books  on  the  subject 
which  include  Minieka  (1979),  Kennington  and  Helgason  (1980)  and 
Jensen  and  3ames  (1980). 

Along  with  the  computational  advances,  network  models  have 
been  applied  to  a  wide  range  of  problem  situations.  In  particular, 
several  water  resource  applications  are  reported  by  the  Texas 
Department  of  Water  Resources  (Texas  Water  Development  Board 
(1974a,  1975))  and  Jensen  et  al.  (1974).  This  report  deals  with 
the  application  of  a  new  network  model  formulation  as  applied  to  a 
water  resources  system. 

Optimal  operation  of  a  system  of  interconnected  water 
reservoirs  i3  an  important  problem  in  water  resources  management. 
The  limited  water  resources  available  coupled  with  the  diverse, 
often  competitive,  projected  demands  on  these  resources  appear  to 
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place  potentially  unacceptable  limits  upon  the  achievement  of 
economic,  social  and  environmental  goals. 

The  operation  of  a  muliperiod  multireservoir  system 
requires  that  the  system  controller  make  decisions  regarding  the 
storage  or  release  of  water  for  each  of  the  reservoirs  on  a 
periodic  basi3.  This  period  may  be  daily,  weekly,  monthly,  etc. 
His  decisions  may  be  based  upon  the  amount  of  water  available  to 
him  in  each  of  the  reservoirs,  the  types  of  demands  for  water  from 
the  various  users  and  upon  his  anticipation  of  the  future 
availability  of  water.  Each  period  lends  itself  to  a  network 
representation  similar  to  the  one-period,  three- reservoir  system 
shown  in  Figure  1 -1 . 

The  amount  of  water  available  for  distribution  i3  a 
function  of  the  amount  of  water  stored  from  the  previous  period, 
the  amount  of  inflow  from  runoff  or  from  upstream  releases  within 
the  period  and  any  purchases  from  outside  sources.  Most  work  to 
date  has  treated  the  amount  of  water  from  runoff  and  other  flows 
into  or  out  of  the  system  deterministically.  That  is,  all  data  and 
parameters  of  the  models  are  assumed  to  be  known  with  certainty. 
Thus  the  models  represent  a  decision  problem  in  which  the  decision 
maker  is  faced  with  a  great  deal  of  complexity  but  no  uncertainty. 
The  complexity  makes  the  decision  problem  difficult  in  itself.  If 
there  is  uncertainty  in  the  real  situation  it  is  often  simply 
ignored  by  the  model. 

The  multiperiod  deterministic  model  assumes  that  the  system 


Figure  1-1 

Three  Reservoir  Problem 
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controller  knows  what  the  future  availability  of  water  will  be  with 
certainty.  He  also  assumes  that  the  demands  placed  upon  the  system 
ar9  known.  Vi th  this  information  available  to  him  he  can  then 
determine  an  optimal  3et  of  decisions  for  the  entire  time  horizon. 

The  neglect  of  uncertainty  results  in  unrealistic  solutions 
where  a  major  aspect  of  the  decision  process  relates  to  dealing 
with  uncertainty.  The  water  resources  problem  obviously  is  dynamic 
in  nature  in  that  decisions  must  be  made  sequentially  over  time. 
Uncertainty  play3  a  significant  role  in  the  decision  process  due  to 
the  unpredictabiltiy  of  nature  in  its  supply  of  surface  waters  to 
the  system  and  also  to  the  incomplete  predictability  of  the  actions 
of  man  in  his  use  of  the  available  resource.  It  is  clear  that  the 
controller  of  the  system  must  exhibit  caution  in  setting  reservoir 
levels  and  river  releases  so  that  unlikely  but  possible  natural 
events  do  not  cause  the  system  to  fail  in  its  functions  of 
providing  a  reliable  water  supply  and  protection  against  floods. 

A  deterministic  model  does  not  exhibit  caution  in  a 
dynamic,  multiperiod  model.  Since  all  data  is  assumed  certain,  the 
future  in  all  required  detail  is  known.  An  optimum  solution  can  be 
determined  which  provides  maximum  benefit  at  minimum  cost.  A 
historical  sequence  of  water  runoff  and  demand  data  might  be  used 
to  give  the  model  "realistic"  data.  The  one  aspect  of  the  solution 
procedure  which  is  not  realistic  for  the  multiperiod  model  is  that 
the  deterministic  solution  algorithm  has  the  ability  to  look  ahead 
in  time  and  prepare  for  the  events  which  are  to  occur.  Thus, 
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decisions  obtained  through  the  models  do  not  tend  to  resemble  the 
real  decision  process.  The  decision  maker,  after  all,  does  not 
have  this  "look  ahead"  capability  of  the  algorithm.  He  must  make 
decisions  in  the  face  of  uncertainty  and  revise  them  as  time  goes 
on  and  as  uncertainties  become  resolved.  The  incorporation  of 
uncertainty  into  a  water  resources  model  is  a  part  of  this 
research. 

Another  problem  with  most  existing  models  is  that 
traditional  planning  methodology  has  generally  been  directed  toward 
the  analysis  of  projects  individually  in  an  effort  to  match 
reservoir  operation  with  anticipated  requirements.  When  the 
interaction  of  individual  reservoirs  became  more  pronounced  and 
could  not  be  ignored,  operating  criteria  were  often  still  selected 
on  the  basis  of  these  single-project  analyses  through  coordinated 
single- reservoir  simulation  studies.  It  is  noted  however,  that  in  a 
serially  connected  system  of  reservoirs,  the  value  of  water  stored 
in  a  particular  reservoir  is  affected  by  the  amount  of  water  stored 
in  other  reservoirs.  This  relationship  has  been  neglected  in  the 
past  due  to  the  fact  that  this  interaction  between  reservoirs 
suggests  a  nonseparabl9  benefit  function  as  a  function  of  all  the 
reservoirs  versus  a  separate  benefit  function  for  each.  This  means 
that  the  total  benefit  of  the  system  cannot  be  measured  simply  by 
summing  the  individual  reservoir  benefits  as  a  function  of  their 
contents.  One  reason  for  this  neglect  lies  in  the  difficulty  in 
determining  with  any  confidence  just  what  this  joint  function 
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3hould  be.  The  benefit  function  that  will  be  used  herein  is  a 
nonseparabls  quadratic  function  which  measures  or  reflects  the 
current  and  future  value  of  water  stored  by  the  system.  Thus, 
besides  individual  reservoir  benefits,  the  interactive  or  .joint 
reservoir  benefits  will  also  be  evaluated..  This  idea  combined 
with  the  multiperiod  decision  process  is  used  in  a  dynamic 
programming  approach  to  successively  generate  these  benefit 
functions. 

Thi3  report  describes  a  method  to  overcome  the  deficiencies 
of  the  deterministic  solutions  while  including  the  provision  for 
evaluating  interactive  reservoirs,  network  models  are  still  used 
but  the  model  is  changed  in  such  a  way  as  to  exhibit 
characteristics  of  the  true  decision  process.  Full  advantage  i3 
taken  of  the  network  structure  of  the  problem  by  utilizing 
extensively  the  computational  techniques  that  have  been  so 
successful  for  deterministic  models.  Embedding  this  network  model 
in  a  dynamic  programming  solution  approach  which  begins  at  some 
specified  and  finite  future  date  and  works  backward  in  time  to  the 
present  provides  the  necessary  data  to  allow  the  derivation  of 
successive  benefit  functions  which  reflect  the  future  value  of  a 
given  configuration  of  reservoir  contents. 

Chapter  3  describes  the  deterministic  network  models  and 
provides  a  brief  survey  of  the  computational  techniques  used  to 
solve  them.  Chapter  3  also  provides  the  notational  basis  for  the 
remaining  chapters  of  the  report. 
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In  Chapter  4,  the  dynamic  programming  algorithm  for  solving 
the  multiperiod  multireeervoir  stochastic  problem  will  be 
presented.  Chapter  5  includes  the  nonlinear  network  solution 
methodology.  These  network  solutions  embedded  in  a  dynamic 
programming  methodology  provide  the  basis  for  deriving  a  functional 
representation  of  the  future  value  of  water  in  the  face  of 
uncertainty.  In  Chapter  6  the  statistical  aspects  of  the  problem 
will  be  addressed. 

Chapter  7  includes  some  example  applications  which  are 
supported  by  data  contained  in  the  Appendix. 

1.2  Primary  Contributions  of  This  Research 

One  contribution  of  this  research  is  the  development  of  a 
generalized  network  model  which  is  capable  of  solving  network  type 
problems  where  some  of  the  arcs  have  nonlinear  quadratic  cost 
functions.  These  functions  are  allowed  to  be  nonseparable  and  the 
model  is  solved  without  reverting  to  piecewise  approximations  for 
the  arc  costs.  The  only  restriction  is  that  the  overall  objective 
function  which  is  a  combination  of  several  linear  terms  along  with 
some  quadratic  terms  be  a  CONVEX  cost  function;  or  a  CONCAVE 
benefit  function  in  this  case.  It  is  noted  that  there  exists 
several  other  techniques  which  could  be  used  for  this  class  of 
problems.  Some  of  these  are  listed  here: 

1 .  Prank -Wolfe  Method 


2.  Convex  Simplex  Method 
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3.  Method  of  Feasible  Directions 

4.  Gradient  Projection  Methods 

5.  Quadratic  Programming  Algorithm 

9.  Reduced  Gradient  Method 

7.  Newtons  Method 

9.  Steepest  Descent  Method 

9.  Variable  Metric  Methods 

-  Davidon  Fletcher  Powell 

-  3FGS  ( BROYDEN , FLETCHER , GOLD? ARB , 3HANN0 ) 

Most  of  these  methods  could  be  used  either  directly  or  in  a 
specialized  manner  for  this  class  of  problems.  The  choice  of 
introducing  network  theory  as  still  another  way  to  solve  problems 
of  this  nature  both  expands  and  enhances  the  power  of  network 
theory  as  well  as  providing  an  alternative  to  the  above  3uggessted 
methods. 

The  second  and  primary  contribution  of  this  research 
involves  the  integration  of  this  network  solution  technique  into  a 
much  larger  dynamic  programming  model.  This  larger  model  is  used 
to  solve  multireservoir  multiperiod  water  resources  problems  in  the 
presence  of  uncertain  inflows.  Its  uniqueness  li9S  not  only  in  its 
use  of  a  new  network  subproblem,  but  also  in  the  manner  in  which  it 
utilizes  future  stochastic  runoff  information  and  recursively 
generates  benefit  functions  which  represent  the  net  current  and 
future  expected  benefit  to  the  system  as  a  function  of  the  observed 
current  water  levels.  This  functional  representation  allows  any 
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level  to  be  evaluated  for  determining  the  current  decisions.  Thi3 

i 

form  represents  a  continuous  spectrum  for  current  levels  in  lieu  of 
the  more  common  discretized  levels  in  standard  dynamic  programming 
algorithms.  This  means  that  the  current  levels  are  not  required  to 
be  equal  to  or  rounded  to  the  nearest  discrete  level  which  would 
induce  error  into  the  results.  Although  a  discretization  scheme  is 
used  to  derive  the  benefit  functions,  it  i3  done  30  only  to  gain  a 
representation  of  the  true  return  "functionally".  Once  this 
function  is  available,  any  and  all  reservoir  levels  can  be 
evaluated . 

One  main  advantage  of  this  functional  approach  is  that  it 
greatly  relieves  the  dimensionality  problems  associated  with 
discrete  representation  of  large  dynamic  programming  problems. 

In  the  end,  the  result  is  a  realistic  and  usable  water 
resources  model  since  it  does  account  for  the  uncertainties  of  the 
future.  It  can  handle  larger  models  due  to  the  functional 
approach,  and  perhaps  most  importantly,  the  entire  model  has  been 
implemented  into  a  workable  computer  program  where  it  is  readily 
available  to  such  users  as  the  Texas  Water  Development  Board. 

« 
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CHAPTER  II 


2.  LITERATURE  REVIEW 


2.1  General 

Host  of  the  research  involving  analytical  modeling  of 
multireservoir  water  systems  has  occured  during  the  last  15  years. 
There  has  been  a  great  deal  of  variation  in  the  mathematical 
techniques  employed.  Roefs  (1963),  3uras  (1972),  and  Hall  and 
Dracup  (1970)  discuss  the  mathematical  techniques  used  and  the 
variations  of  the  problem  for  which  each  technique  is  most  suited. 
Butcher  ( 1 977 )  indicates  that  a  mathematical  tool  useful  for  one 
water  resource  problem  may  not  be  suitable  for  other  seemingly 
similar  problems.  Multireservoir  models  can  be  roughly  classified 
into  three  categories  depending  upon  their  emphasis  and  scope. 

1.  Design  Models.  These  types  of  models  make  decisions 
concerning  the  construction  of  the  reservoir  system.  They  are 
sometimes  called  capacity  expansion  models.  Decisions  are  made 
concerning  the  size,  location,  and  time  of  construction  of 
reservoirs  and  canals  in  addition  to  determining  water  allocation. 

2.  Water-use  Models.  In  these  models  the  reservoirs  are 
considered  to  be  multipurpose;  that  is,  several  possible  uses  of 


water  are  available  at  each  reservoir. 


Decisions  are  made 


concerning  such  things  as  the  timing  and  extent  of  irrigation  for 
various  crops.  These  models  most  often  concern  a  single  reservoir. 

3-  Time  Planning  Models.  The  main  objective  with  these 
models  is  to  determine  the  use  and  storage  of  water  in  several 
interconnected  reservoirs  in  3uch  a  way  as  to  be  prepared  for 
future  shortages  or  excesses.  Thi3  research  concerns  itself  with 
this  kind  of  model. 

Typically,  literature  regarding  water  reservoir  operations 
is  characterized  by  four  primary  factors.  These  being: 

1.  System  -  single  versus  multireservoir 

2.  Operation  -  single  versus  multiperiod 

3.  Inflows  -  Deterministic  versus  stochastic 

4.  Return  or  objective  function  -  Linear  or  separable 

nonlinear  versus  nonlinear  nonseparable. 

The  mathematical  models  employed  to  model  reservoir 
problems  have  included  the  following: 

1 .  Linear  programming 

2.  Dynamic  programming  -  both  deterministic  and  stochastic 

3.  Chance-constrained  linear  programming 

4.  Decomposition  approaches 

5.  Simulation 

5.  Markov  chains 

7.  Networks 

S.  Nonlinear  programming 
In  most  cases  combinations  of  the  above  were  used. 
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Many  authors  have  used  the  above  techniques  to  develop 
models  for  single  reservoir  or  single  period  systems  applying  both 
deterministic  and  stochastic  inflows.  These  works  are  presented 
briefly  in  section  2.2.  Since  this  research  is  concerned  with 
multireservoir  multiperiod  systems,  a  finer  breakdown  of  models  as 
they  apply  to  the  multireservoir  multiperiod  systems  is  discussed 
in  sections  2.3,  2.1.  and  2.5. 

2.2  Single  Reservoir  (single  and  multi  period)  or  Multireservoir 
Single  Period  Models: 

Network  models  include  the  techniques  used  by  the  Texas 
Water  Development  Board  (1974a,  1974b)  and  those  mentioned  below. 
3haumik  (1973)  presents  an  optimum  operating  policy  of  a  water 
distribution  system  with  losses  (gains  less  than  one).  Concern 
here  was  with  the  economic  effect  of  seepage  and  evaporation  of 
water  from  canals  and  reservoirs  in  a  water  distribution  sytem  with 
reference  to  the  Texas  water  plan.  Weirs  and  Beard  (1971)  and 
Svenson  and  Mosely  ('975)  present  more  detailed  discussions  of  the 
Texas  Water  Development  Board's  work  in  time  planning  and  design 
models . 

Linear  programming  has  been  applied  to  all  types  of 
reservoir  problems.  ReVelle  and  Gundelach  (1975)  introduced  a  new 
version  of  the  linear  decision  rule  in  '975.  This  new  form  permits 
the  minimization  of  the  3um  of  the  variances  of  releases,  a 
performance  objective  not  previously  subject  to  the  control  of  the 
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designer.  The  minimization  of  release  variances  of  a  single 
reservoir  allows  further  diminishment  of  losses  associated  with 
deviations  from  target  releases.  They  concluded  that  thi3  new 
formulation,  while  experiencing  definite  advantages  with  regards  to 
the  minimum  release  level,  had  a  disadvantage  in  that  it  required 
larger  reservoir  capacities.  In  spite  of  this  result,  this  new 
linear  decision  rule  makes  it  possible  to  attain  release  levels 
that  otherwise  might  be  regarded  as  infeasible,  Sundelach  and 
SeVelle  (1975)  then  use  this  new  decision  rule  and  develop  a 
chance-constrained  model  which  seeks  the  smallest  reservoir 
satisfying  certain  conditions  on  storage,  release  and  freeboard. 

These  linear  programming  methods  are  considerably  slower 
than  the  linear  out-of-kilter  algorithm,  but  they  have  more 
flexibility.  For  instance,  monthly  evaporation  can  be  included  in 
the  model  as  a  percentage  of  reservoir  storage.  For  a  multiperiod 
model,  linear  programming  has  the  shortcommings  of  the  necessity 
for  perfect  information,  the  large  size  of  the  problem,  failure  to 
make  use  of  the  final  reservoir  storage,  wasted  computation  and  of 
course  the  restriction  of  linearity.  As  stated  by  Buras  (1972), 
"linear  programming  yields  only  point  solutions  in  the  policy 
space,  no  matter  how  many  dimensions  the  space  has.  Most 
situations  in  which  the  state  of  the  system  changes  (in  time  or  in 
space)  and  in  which  decisions  have  to  be  taken  successively  are 
clearly  outside  the  grasp  of  linear  programming" .  A  point  solution 
means  the  set  of  optimal  values  for  each  of  the  variables,  given 


fixed  values  for  each  of  the  parameters  of  the  3ystem.  This 
differs  from  a  functional  type  of  solution  where  the  solution 
variables  are  given  as  a  function  of  another  variable. 

Dynamic  programming  is  the  most  theoretically  appealing 
approach  to  multiperiod  reservoir  models  of  all  types  since  these 
problems  involve  sequential  decision-making  processes.  Also,  the 
outcome  of  each  decision  (or  set  of  decisions  in  a  time  period) 
appears  as  a  function  rather  than  as  a  point  solution.  That  is, 
the  optimal  decision  to  be  made  i3  determined  for  any  state  of  the 
system.  This  allows  3uboptimal  policies  to  be  examined,  a 
desirable  feature  due  to  the  inherent  uncertainty  in  multireservoir 
problems.  This  feature  makes  dynamic  programming  especially  useful 
for  real  time  system  operation.  The  limitation  on  the  usefulness 
of  dynamic  programming  is  the  so  called  "curse  of  dimensionality". 
Each  reservoir  gives  rise  to  a  new  3tate  variable  (usually  the 
final  reservoir  storage  level).  If  there  are  R  reservoirs  and  each 
reservoir  has  K  possible  levels,  there  are  K®  possible  state 
combinations  per  time  period.  For  this  reason  most  of  the  dynamic 
programming  models  have  been  for  systems  with  either  one  or  two 
reservoirs. 

Ruras  (1972)  fixes  four  or  five  as  the  maximum  number  of 
reservoirs  that  can  be  handled  computationally  by  dynamic 
programming.  Dynamic  programming  al30  fails  to  take  into  account 
the  stochastic  nature  of  the  multireservoir  problems.  Thi3  can  be 
rectified  by  using  stochastic  dynamic  programming  as  was  done  by 
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Driscoll  (1974-). 

3uras  (1972)  presents  a  dynamic  programming  formulation  for 
a  design  type  model  and  for  a  water-use  model.  Young  (1967) 
combines  dynamic  programming  with  a  simulation  approach. 

Butcher  (1973)  defines  stochastic  dynamic  programming  to  be 
those  formulations  of  dynamic  programming  in  which  the  value  of  one 
of  the  state  variables  is  related  in  a  probabilistic  way  to  the 
value  of  that  3ame  variable  in  adjacent  time  periods.  In  terms  of 
multireservoir  problems,  stochastic  dynamic  programming  allows  the 
rainfall  and  demand  (or  net  demand,  ie.  demand  minus  inflows)  at  a 
reservoir  in  one  time  period  to  be  dependent  probabilistically  on 
the  net  demand  at  that  reservoir  in  the  previous  time  period.  The 
optimal  policy  developed  is  that  which  minimizes  expected  costs  for 
the  system.  As  with  traditional  dynamic  programming,  the  optimal 
policy  is  in  a  form  that  can  readily  be  used  for  real  time  system 
operation.  However,  the  dimensionality  problem  is  compounded  due 
to  the  additional  state  variables  which  are  the  net  demands  at  the 
various  reservoirs  in  the  previous  period. 

Due  to  the  dimensionality  problem,  stochastic  dynamic 
programming  models  have  been  applied  mainly  to  water-use  models, 
rather  than  other  models  which  tend  to  have  more  than  one 
reservoir. 

Butcher  (1971)  presents  such  a  model  for  one  multipurpose 
reservoir.  Loucks  (1969)  presents  three  stochastic  dynamic 
programming  models  that  he  used  to  define  operating  policies  for 
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several  of  the  ?inger  lakes  in  Jfew  York.  These  are  also  one 
reservoir  models.  Instead  of  economic  objectives,  Loucks  minimizes 
the  sum  of  squares  of  the  departures  of  releases  from  a  set  of 
target  releases  specified  by  the  state. 

Chance  constrained  linear  programming  is  another  tool  that 
has  been  applied  to  multireservoir  problems  in  an  attempt  to 
account  for  stochastic  variation.  Like  deterministic  linear 
programming,  thi3  method  is  more  suitable  for  design  or  water-use 
models  than  it  is  for  time  planning  models.  The  fact  that  point 
solutions  are  found,  causes  chance  constrained  linear  programming 
to  be  less  applicable  to  real-time  system  operation  over  a  long 
time  span.  Loucks  (1969)  proposes  a  one-reservoir  water-use  model. 
ReVelle  at  al.  (1976)  considered  the  use  of  linear  decision  rules 
and  the  development  of  a  stochastic  model  in  1969.  In  1975,  Loucks 
and  Dorfman  ('969)  compared  several  chance  constrained  linear 
decision  models  for  reservoir  planning  and  operation.  Their  basic 
conclusion  was  that  while  all  results  of  the  four  decision  rules 
considered  were  within  the  constraints  of  the  problem,  all  tend  to 
yield  overly  conservative  results.  They  state  that  linear  decision 
rules  permit  the  use  of  linear  programming  methods  for  solving  what 
would  otherwise  be  a  very  messy  nonlinear  stochastic  optimization 
problem.  This  is  indeed  a  mathematical  advantage,  but  at  the  same 
time,  these  linear  decision  rules  reduce  considerably  the  number  of 
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possible  operating  policies  that  can  be  considered.  Hence,  the 
rule  itself  is  a  constraint  imposed  for  mathematical  reasons. 

Klemes  (1977)  and  Doran  (1975)  discuss  the  problem  of 
selecting  discrete  reservoir  levels  and  how  the  "curse  of 
dimensionality"  can  be  overcome  by  using  a  method  called  the 
divided  interval  technique.  This  technique  differs  from  the 
traditional  discretization  method  primarily  in  the  precision  with 
which  the  two  boundary  states  are  represented.  The  traditional 
method  developed  by  Moran  (1954)  tends  to  over  estimate  the 
probabilities  of  emptiness  and  fullness,  thereby  underestimating 
the  intermediate  levels.  Klemes  and  Doran  show  that  for  equivalent 
results,  5-10  discretizations  using  the  divided  interval  technique 
corresponds  to  approximately  30  intervals  using  the  traditional 
method. 

The  remaining  mathematical  methods  have  been  less  used  and 
are  not  easily  classified.  Parikh  (1966)  and  (1967)  presents  a 
linear  decomposition  method  designed  for  use  in  a  northern 
California  system.  Young  (1967)  combines  dynamic  programming  with 
Monte-Carlo  simulation  of  stochastic  inflows.  Su  and  Deininger 
(1971)  present  a  Markov-chain  approach  for  serially  connected 
reservoirs.  He  solves  the  Markov  system  by  a  method  of  successive 
approximations  rather  than  by  dynamic  programming. 

2.3  Deterministic  Inflows  and  Linear  or  Separable  Nonlinear 


Objective  Function: 
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Within  this  category  of  multiperiod  multireser/oir  systems 
many  techniques  have  been  used.  Several  authors  used  a  linear 
programming  approach.  Drobny  (1971)  was  concerned  with  water 
quality  and  quantity  problems.  Salcedo  (1972)  dealt  with  a 
water-use  model.  Mannos  (1955)  was  concerned  with  the  efficiency 
of  operation  of  a  system  of  iams  and  Me.jia  et  al.  (1974)  evaluated 
multireservoir  operating  rules  using  linear  programming. 

Schweig  and  Cole  (1968)  used  dynamic  programming  to  deal 
with  random  inflows  having  first  order  serial  correlation.  The 
distribution  of  these  random  inflows  was  approximated  by  discrete 
probability  space.  This  correlation  was  simplified  by  classifying 
the  inflow  data  according  to  whether  an  item  was  preceded  by  an 
inflow  higher  or  lower  than  mean  for  the  antecedent  month.  Rood  at 
al.  (1975)  presents  a  dynamic  programming  model  for  the 
time-planning  variety  that  i3  especially  designed  for  serially 
linked  reservoirs  thereby  reducing  the  state  space  dimensionality 
problem. 

Fults  and  Hancock  (1972)  use  state  incremental  dynamic 
programming  to  find  the  optimal  operating  policy  for  a  four 
reservoir  system.  The  objective  is  to  maximize  power  generation 
while  satisfying  firm  water  contracts,  enhancing  environmental 
aspects  and  providing  flood  control. 


Heidari  et  al.  (1971)  and  Meredith  (197b)  use  a  technique 
called  discrete  differential  dynamic  programming  (DDDP).  This  is  an 
iterative  method  that  eases  the  state  dimensionality  problem  by 


starting  with  a  trial  trajectory  satisfying  a  specific  set  of 
initial  and  final  conditions  and  applies  Bellman's  recursive 
equation  in  the  neighborhood  of  this  trajectory.  At  the  end  of 
each  iteration  a  logically  improved  trajectory  is  obtained  and  used 
as  the  trial  trajectory  in  the  next  step. 

Becker  and  Yen  (1974)  use  a  combination  of  linear  and 
dynamic  programming  for  the  optimization  of  real  time  operations  of 
a  multi reservoir  system  and  Hirsch  et  al.  (1977)  combine  linear 
programming  with  simulation  techniques. 

Prekopa  et  al.  (1968)  address  serially  linked  reservoir 
design  by  attempting  to  meet  all  demands  with  a  given  high 
probability.  Their  objective  is  to  minimize  the  3um  of  the 
building  costs  and  penalties  incurred  for  unsatisfied  demand. 
Their  method  of  solution  uses  a  sequential  constrained  minimization 
techni^ -e  '5UMT)  with  a  logarithmic  penalty  function.  A 
disadvantage  of  this  approach  is  that  under  certain  circumstances 
not  all  demands  can  be  met  with  the  desired  probability. 

Jensen  et  al.  1T974)  represent  a  multireservoir  multiperiod 
system  using  networks.  Here,  the  network  for  each  period  stays  the 
3ame  with  inclusion  of  a  3et  of  storage  arcs  to  join  the  networks 
from  one  period  to  the  next.  Kerr  C 1 972 )  combines  linear 
programming  and  the  out-of-kilter  algorithm  and  applies  these  to 
the  Saskatchewan-Helson  river  basin  in  Canada.  He  compares 
multireservoir  analysis  techniques  by  considering  53  possible 
future  storage  reservoirs  and  22  diversion  possibilities. 
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Hirsch,  Cohon  and  ReVelle  (1977)  developed  a  hypothetical 
design  for  the  3izing  of  three  reservoirs  in  parallel.  They 
recognized  that  there  are  benefits  due  to  the  .joint  operation  of  a 
system  of  reservoirs  in  excess  of  the  benefits  from  optimal 
individual  operation.  Basically,  they  concluded  that  within 
reasonable  limits  any  combination  of  three  reservoirs  whose 
capacities  3um  to  the  same  total  capacity  has  nearly  the  same 
maximum  system  yield.  Their  objective  wa3  to  meet  all  demands, 
which  were  deterministic,  while  not  accounting  for  spillage  or 
evaporation.  The  method  of  solution  involved  simulation  of  five 
years  of  actual  data  and  a  linear  programming  optimization  model. 

Windsor  and  Chow  (1972)  present  a  mixed  linear  programming 
model  with  integer  variables  that  is  a  combination  design  and 
water-use  model.  The  linear  programming  method  appears  to  be  more 
suitable  for  models  of  these  types  where  the  number  of  time  periods 
can  be  kept  to  a  minimum.  They  consider  their  model  a  practical 
one  computationally ,  but  admit  it3  weakness  in  not  considering  the 
stochastic  nature  of  the  problem,  tfeier  and  Beightler  (1967)  use  a 
decomposition  method  for  branching  multi-stage  water  resource 
systems. 


In  the  area  of  separable  nonlinear  objective  functions,  Lee 
and  Waziruddin  (1970)  use  two  approachs,  the  gradient  projection 
and  conjugate  gradient  methods.  They  consider  the  profit  accrued 
from  irrigation  and  the  benefit  received  from  recreation  to  be 
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Roefs  and  Bodin  (1970)  use  separable 


functions. 


programming  with  Dantzig-'Wolfe  decomposition. 

2.4.  Deterministic  Inflows,  Non-Separable  Objective  Function 

In  this  category,  Gagnon  et  al.  (1974)  use  a  generalized 
reduced  gradient  approach  to  a  very  large  hydroelectric  system. 
Lasdon  (1976)  and  TVA  (1974)  also  use  generalized  reduced  gradient 
methods  as  applied  to  water  resource  systems.  Trott  and  Yeh  (1973) 
use  a  3tepwise  state  variable  incremental  dynamic  programming 
approach  and  TVA  (I974a)  used  a  dynamic  programming  successive 
approximation  approach. 

Lui  and  Tedrow  (1973)  use  dynamic  programming  and  a  multi 
variable  pattern  search.  This  multi  variable  polynomial  objective 
function  represents  the  current  and  future  economic  losses  to  the 
system.  Thi3  function  is  determined  by  regression  analysis  where 
the  states  or  reservoir  levels  are  the  independent  variables  and 
the  functional  return  is  the  dependent  variable.  This  method  of 
representing  future  economic  losses  functionally  tends  to  eliminate 
the  problem  of  dimensionality.  They  use  a  random  sampling 
technique  to  assure  unbiased  and  efficient  selection  of  initial 
state  variable  level  combinations.  In  Rosenthal  (1977),  a 
multiperiod  multireservoir  release  scheduling  for  maximum 
hydropower  benefit  was  formulated  by  the  Tennessee  Valley  Authority 
as  an  optimization  model  with  a  nonseparable  nonlinear  objective 
function  and  linear  network  constraints.  Rosenthal  presents  a  new 
solution  technique  based  on  reduced  gradient  techniques  and  on 
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primal  linear  network  flows.  An  unusual  feature  of  the  algorithm 
is  an  integer  programming  subproblem  whose  exact  solution 

determines  the  search  directions.  Test  problems  were  run  on  a  six 
reservoir  TVA  system.  His  network  is  somewhat  unique  and 

constrained  in  that  the  system  is  required  to  be  an  arborescence . 

An  arborescence  i3  a  tree  with  the  property  that  no  two  arcs  are 
directed  away  from  the  same  node,  and  a  tree  is  a  connected 

loopless  network.  Thus,  in  thi3  case,  no  provisions  are  made  for 
pipeing  or  channeling  water  to  other  locations.  All  water  flows 
downstream  to  the  next  reservoir  in  series. 

2.5  Stochastic  Inflows  and  Linear  or  Separable  Nonlinear  Objective 
Function 

Consideration  of  the  stochastic  nature  of  inflows  to 
reservoir  systems  for  multireservoir  multiperiod  systems  has  just 
recently  begun  to  attract  attention.  Sobel  (1975)  analyzes  the 
structure  of  optimal  policies  for  several  discrete  time  control 
models  of  reservoir  storage  using  dynamic  programming.  Most  of  the 
models  considered  are  stochastic  3nd  are  prompted  by  operating 
problems  of  regulating  the  amounts  of  water  discharged  from 
reservoirs.  He  develops  an  analogy  between  models  of  multiple 
reservoir  systems  and  of  multi-item  inventory  models.  Pinter 
('976)  uses  a  stochastic  dynamic  programming  method.  Driscoll 
(1974)  uses  a  stochastic  dynamic  programming  approach  to  a 
multireser/oir  multiperiod  problem  that  uses  a  revised  nonlinear 
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out-of-kilter  algorithm  developed  by  Jensen  and  Reeder  (1974)  at 
each  stage  storing  the  results  in  a  benefit  matrix.  3y  assuming  a 
cyclical  pattern  he  repeatedly  cycles  through  the  periods  until 
convergence  of  the  benefit  function  is  attained.  Thi3  revised 
out-of-kilter  algorithm  allows  convex  functions  and  a  limited  type 
of  nonseparable  cost  functions.  Thi3  method  allowed  him  to 
consider  systems  with  five  reservoirs  without  too  much  difficulty. 

Thu  (1980)  developed  a  method  to  deal  with  stochastic 
situations  with  recourse.  His  work  involved  a  two  3tage  decision 
process  whereby  an  initial  decision  was  mads  based  on  expected 
demands  and  then  as  actual  demands  became  known,  a  second  decision 
(the  recourse)  was  made  to  satisfy  all  demands.  His  objective  wa3 
to  minimize  the  sum  of  the  costs  from  the  first  decision  and  the 
expected  penalty  cost3  as  required  by  the  recourse  actions. 

Roefs  0  963)  presents  a  stochastic  dynamic  programming 
formulation  for  one  and  two  reservoir  systems. 

Turgeon  (i960)  uses  two  mathematical  manipulation 
techniques  to  solve  problems  too  large  for  dynamic  programming. 

1.  The  one  at  a  time  method,  which  breaks  up  the 
multi-variate  problem  into  a  series  of  one  state  variable 
subprobleras,  and 

2.  Aggregation/decomposition  method  which  breaks  up  the  N 
state  variable  problem  into  H  subproblems  of  two  state 
variables. 

Both  of  these  methods  are  then  solved  using  dynamic  programming. 


He  applies  these  methods  to  a  six  reservoir  system  where  the 
inflows  are  assumed  normal  with  a  mean  and  variance  corresponding 
to  those  of  the  Quebec  river  historical  data. 

Joeres  et  al.  (1971)  considered  chance  constrained  linear 
programming  in  conjunction  with  linear  programming  techniques  in 
deriving  operating  rules  for  joint  operation  of  raw  water  sources. 
Curry  and  Helm  (1972)  present  a  chance  constrained  model  for  a 
3ingls  multipurpose  reservoir  and  then  Curry  et  al.  (1972)  extend 
this  to  a  system  of  linked  multipurpose  reservoirs.  They  allow  the 
unregulated  inflows  into  each  reservoir  at  each  time  period  to  be 
stochastic  with  a  known  probability  distribution.  There  is 
independence  between  the  reservoirs  and  for  the  same  reservoirs  in 
different  time  periods.  They  show  how  their  formulation  reduces  to 
a  deterministic  linear  programming  model. 

Helm,  Curry  and  Hasan  (1972)  present  a  design  model  for  a 
system  of  reservoirs.  This  formulation  employs  a  mixed  continuous 
and  integer  linear  programming  form  that  is  solved  by  Benders 
decomposition  method. 

Sigvaldason  (1975)  uses  simulation  and  the  out-of-kilter 
algorithm  and  applies  his  model  to  the  Trent  River  system  in 
Ontario,  Canada.  He  divides  each  reservoir  into  five  storage  zones 
and  applies  penalty  coefficients  for  any  deviations  from  ideal 
conditions  as  applied  to  these  zones.  3odin  and  Roefs  (1971 )  use 
separable  programming  and  the  Oantzig-Wolfe  decomposition  method. 

Houck  and  Cohon  (1973)  utilized  a  sequentially  explicitly 
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3tochastic  linear  programming  model  (SELSP)  to  determine  a  design 
and  management  policy  for  a  two  dam  system.  Their  process  was  to 
sequentially  solve  policy  and  design  linear  programming  models 
which  are  formed  by  specifying  minimal  portions  of  the  nonlinear 
program.  The  major  wealcnesses  of  the  model  are  high  data 
requirements  and  computational  burden.  They  propose  a  method  of 
mitigating  both  deficiencies  while  explicitly  retaining  the 
interaction  of  the  reservoir  system.  The  system  coordinated 
performance  individual  operation  (SCORPIO)  method  provides  the 
necessary  information  to  evaluate  the  interaction  among  facilities 
in  a  multi reservoir  system.  It  i3  a  way  to  utilize  the  available 
data  efficiently  and  its  use  makes  the  use  of  SELSP  models 
practical.  Basically,  SCORPIO  solves  the  individual  reservoir 
design  and  operating  problems,  and  then  system-wide  performance 
characteristics  are  obtained  by  using  expected  values  and 
correlations  between  streamflows  at  all  of  the  3ites. 

Takeuchi  and  Moreau  (l97i)  use  a  combination  of  linear 
programming,  dynamic  programming  and  regression  analysis.  The 
monthly  operating  decisions  are  given  by  solution  of  a  piecewise 
linear  program,  the  objective  function  for  which  consists  of  two 
part3.  One,  the  immediate  economic  losses  within  the  month,  and 
two,  the  expected  present  value  of  future  losses  as  a  function  of 
end  of  month  storage  levels.  These  expected  losses  are  determined 
by  imbedding  the  linear  program  in  a  stochastic  dynamic  program. 
Their  loss  functions  were  constructed  in  accordance  with  the 


Si  ■ 


following: 

1 .  Diminishing  marginal  utility  of  water 

2.  Deficits  involving  high  proportions  of  nominal  municipal 
water  use  were  considered  catastrophic 

3.  High  rates  of  deficit  in  low-flow  augmentation  also 
created  serious  damage. 

They  applied  their  model  to  a  five  reservoir  system  in  the  Piedmont 
region  of  North  Carolina. 

2.5  Summary  of  Literature 

In  summary,  none  of  the  literature  reviewed  addressed  the 
full  combination  of  multireservoir,  multiperiod,  stochastic  inflows 
and  nonseparable  objective  functions.  Rosenthal  (1980)  made  this 
same  observation  after  having  researched  over  100  articles. 

The  chapters  to  follow  present  a  new  approach  for  this 
combined  set  of  conditions.  The  methods  employed  primarily  use 
networks,  dynamic  programming  and  regression  analysis.  The  manner 
in  which  these  are  used  encompasses  the  stochastic  nature  of 
inflows  using  a  Monte  Carlo  approach. 


CHAPTER  III 


3.  The  Deterministic  Model 

3.1  The  Network  Model 

This  chapter  is  used  to  describe  the  network  flow  model, 
introduce  the  notation  to  be  used  throughout  the  report,  and  review 
the  solution  approaches  used  to  solve  deterministic  problems.  The 
latter  are  used  extensively  in  the  algorithms  which  solve  the 
stochastic  problem  as  well.  This  chapter  is  taken  from  Chapter  II 
of  Jensen  et  al.  (1980),  co-authored  by  the  author  of  this  report. 

A  network  flow  model  is  simple  in  the  sense  that  it 
requires  very  few  kinds  of  structural  elements  and  parameters  to 
describe  it.  A  complex  model  is  constructed  from  imaginative 
arrangements  of  these  simple  elements.  Figure  3-1  illustrates  the 
basic  structure  of  a  network.  The  network  consists  of  nodes  and 
arcs.  The  nodes  are  represented  by  circles  with  the  inscribed 
number  used  for  identification.  For  general  reference  lower  case 
letters  such  as  i  and  j  are  used  to  refer  to  nodes,  where  the 
letters  symbolize  numeric  identifiers.  Arcs  are  the  directed  line 
segments  going  from  one  node  to  another.  Consider  the  general  arc 
k  (k  refers  to  a  numeric  identifier  assigned  to  the  arc) 
originating  at  node  i  and  terminating  at  node  j.  Frequently  the 
notation  k(i,j)  is  used  in  cases  where  it  is  important  to  emphasize 


27 


V, 


29 


» 


« 


« 


the  identity  of  the  nodes  touching  arc  It. 

The  variable  quantities  associated  with  the  network  are  arc 
flows.  The  symbol  f  represents  the  flow  in  arc  k.  The 
optimization  problem  is  to  determine  the  values  of  f  for  each  arc 
which  minimize  some  criterion  subject  to  certain  constraints.  This 
criterion  may  be  cost,  time,  distance,  etc.  The  criterion  and 
constraints  are  defined  below. 

Associated  with  each  arc  k  are  four  parameters:  lower  bound 
on  flow,  c^;  upper  bound  on  flow,  c^;  marginal  cost,  h^;  and  gain 
a,  .  Parameters  and  variables  associated  with  arcs  are  shown  in 
parentheses  near  the  arc. 

The  values  of  c^  and  c^  provide  simple  bounding  constraints 
for  the  flow  on  arc  k: 


c,  <  f,  <  c, 

-k  —  k  —  k 

The  value  of  h^  indicates  the  marginal  change  in  total  cost  with 

respect  to  f  .  In  linear  problems  is  a  constant  independent  of 

the  value  of  f,  and  the  cost  of  flow  on  arc  k  is: 
k 


A  variation  in  the  form  of  the  arc  cost  function  is  presented  in 
Chapter  5  where  we  will  introduce  a  quadratic  cost  function. 

The  value  of  a.^,  the  arc  gain,  allows  the  flow  to  increase 
or  decrease  as  it  passes  through  the  arc.  For  an  arc  k(i,j)  the 
flow  leaving  node  i  is  f,  .  The  flow  entering  node  j  is  a,  f,  . 
Whether  flow  increases  or  decreases  as  it  passes  through  the  arc 

ieoends  on  the  value  of  a,  .  If  a,  <1 ,  the  flow  decreases.  When 

<  < 


a,  >1  ,  the  flow  will  increase.  If  a,  <0  the  flow  leaving  node  i  with 

tC  < 

value  f,  causes  a  flow  a,  f,  to  leave  node  j  in  the  direction  of 
k  k  < 

node  i  (on  arc  k) .  This  strange  possibility  has  some  applications 

and  is  allowed  by  the  algorithms.  It  ia  only  required  that  a,  WQ 

for  all  arcs,  A  network  in  which  all  gains  are  unity  is  called  a 

pure  network,  while  the  presence  of  one  or  more  nonunity  gains 

results  in  a  generalized  network. 

There  are  also  parameters  associated  with  the  nodes.  Mode 

parameters  and  variables  are  shown  in  brackets  adjacent  to  the 

nodes.  The  most  important  is  node  external  flow,  b^  for  node  i. 

This  parameter  represents  flow  entering  or  leaving  the  network  from 

external  sources  at  node  i.  Use  b  >0  to  imply  flow  entering  node 

i.  If  b^<0,  flow  leaves  the  network  with  magnitude  b^  .  When 

b^O,  no  flow  enters  or  leaves  the  network  at  node  i.  In  the  pure 

network,  the  sum  of  the  flows  entering  the  network  will  equal  the 

sum  of  the  flows  leaving  the  network. 

Additional  node  parameters  are  slack  external  flow,  b  ^ , 

and  slack  external  cost,  h  ..  If  b  .>0,  then  flow  may  be  obtained 

si  si  ' 

from  external  sources  at  node  i  at  a  cost  h  .  per  unit.  The  amount 

si  r 

of  slack  flow,  f  .  is  a  variable  bounded  by  zero  and  b  . .  When 
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b  .  <0,  then  flow  may  be  removed  at  node  i  at  a  cost  of  h  .  per 
unit.  The  amount  of  slack  flow  here  is  bounded  by  zero  and  -b  . 
The  sign  is  used  only  to  indicate  the  direction  of  slack  flow. 
These  slack  external  flows  are  very  useful  modeling  tools. 

One  additional  constraint  type  that  relates  the  flows  in 
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the  network  is  the  requirement  that  flow  be  conserved  at  each  node. 
This  means  that  the  sum  of  the  flows  leaving  a  node  on  the  arcs  of 
the  network  less  the  sum  of  the  flows  entering  the  node  on  network 
arcs  must  equal  the  external  flow  at  the  node. 

The  linear  network  flow  programming  problem  can  be  written 
in  algebraic  form  with  the  definition  of  some  additional  notation 
as  shown  in  Table  3-1  •  With  the  objective  of  minimizing  costs,  the 
problem  is  written  as  follows: 


st 


Model  I 

Minimize  Z  *  hf  +  h  f 

3  3 
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for  i*1  , . .  .n 


c,  <  f.  <  c, 

He  —  k  —  k 

0  <  f  .  <  b  .  for  b  .>0 

—  31  —  SI  31 

0  <  f  .  <  jb  .  I  for  b  .  <  0 

—  31  —  l  31  |  31 

f  .*0  for  b  .  »0 

31  31 
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Notation 

m 

M 

n 


Table  3-1 

Definition  of  Notation 

Definition 

Number  of  arcs  in  the  network 
Set  of  arcs  m*(l ,2,3, . . . ,m) 
Number  of  nodes 
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Set  of  nodes  N*(l ,2,3, . . . ,n) 

Set  of  arcs  that  originate  at  node  i 
Set  of  arcs  that  terminate  at  node  i 
Vector  of  arc  lower  bounds 
Vector  of  arc  capacities 
Vector  of  arc  marginal  costs 
Vector  of  arc  gains 
Vector  of  fixed  external  flows 
Vector  of  arc  flows 
Vector  of  slack  external  flow  bounds 
Vector  of  slack  external  cost3 
Vector  of  slack  external  flows 

For  purposes  of  the  algorithms  to  follow,  two 

transformations  which  allow  a  somewhat  simpler  network  model  are 

now  described.  In  the  algorithms  these  transformations  are 

automatically  performed  by  the  computer  programs.  First  eliminate 

slack  external  flows.  This  i3  done  by  creating  a  new  node  called 

the  slack  node  at  which  conservation  of  flow  is  not  required.  The 

value  of  n  is  increased  by  one  to  account  for  the  slack  node,  and 

this  new  node  is  assigned  the  index  n.  Sow  each  slack  external 

flow  is  replaced  by  an  arc  which  originates  or  terminates  at  the 

3lack  node.  For  each  oositive  value  of  b  ,  at  node  i.  create  an 

si 

arc  from  node  n  to  node  i  with  caoacity  b  ,  and  cost  h  ..  For  each 
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negative  value  of  b  ,  create  an  arc  from  i  to  n  with  capacity 
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l-b  . I  and  cost  h  . .  The  lower  bound  is  set  to  zero  and  the  gain  is 

|  Sll  31 

set  to  one  for  these  arcs.  Now  the  network  has  no  slack  external 
flow  parameters  but  rather  alack  external  flows  are  represented  by 
arc  flows  to  or  from  the  alack  node.  The  arc  set  is  expanded  to 
include  these  new  arcs. 


The  second  transformation  makes  all  arc  lower  bounds  equal 
to  zero.  This  is  illustrated  in  Figure  3-2.  Here,  for  each  arc 
k(i,j)  with  a  nonzero  lower  bound,  c.  ,  adjust  arc  and  node 
parameters  as  follows: 


ck 

■  ^ 

of 

I 

b! 

*  b. 

-  c. 

i 

l 

-k 

'  a 

b  * 

a,  c. 

i 

j 
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The  primed  parameters  are  the  transformed  parameters  which  will  be 

used  by  the  algorithm.  The  effect  of  the  transformation  is  to  make 

all  lower  bounds  zero.  Hence,  they  no  longer  need  be  considered 

explicitly.  To  recover  the  solution  to  the  original  problem,  a 

reverse  transformation  is  required  after  the  problem  is  solved. 

Let  f^  be  the  flow  in  arc  k  obtained  by  the  algorithm  and  f^  be  the 

flow  corresponding  to  the  original  problem.  Then: 

f.*  «f,  *e. 
k  k  -He 

The  cost  of  the  solution  must  also  be  adjusted  accordingly.  Using 
the  same  prime  notation  the  cost  of  arc  k  for  the  original  problem 
is: 


*4 
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With  the  transformed  parameters,  the  minimum  cost 
optimization  problem  now  becomes: 

Model  II 
Minimize  hf 


k€M„ 


a,  f.  *b. 
k  k  i 


i*1 


,n- 


Mote  that  no  conservation  of  flow  constraint  is  written  for  the 
slack  node  as  this  constraint  would  be  redundant. 

Model  II  is  a  bounded  variable  linear  program.  The  matrix 
of  conservation  of  flow  constraints  has  only  two  nonzero  entries 
for  each  arc,  one  equal  to  *1  and  the  other  equal  to  -a^.  When  all 
are  integer  and  all  a^  are  unity  (the  pure  problem)  the  optimum 


flows  will  be  integer.  The  flows  will  not  in  general  be  integer 
for  the  generalized  problem. 
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"5.2  Network  Model  of  the  Multireservoir  System 
3.2.1  The  Multireservoir  System 

One  of  the  principle  ingredients  which  supports  life, 
industry,  agriculture  and  the  environment  in  general  is  water. 
Critical  to  the  continued  preservation  of  life  is  the  intelligent 
and  efficient  use  of  all  available  water  supplies.  This  has 
already  become  a  major  problem  for  many  nations  of  the  world,  and 
fresh  water  supplies  are  becomming  more  and  more  in  demand  as  the 
population  increases  and  as  new  or  improved  industrial  techniques 
require  it.  Water  sources  are  typically  divided  into  ground  water 
and  surface  water,  and  within  the  surface  water  category  they  can 
be  broken  down  into  river  sources,  reservoir  or  lake  sources  and 
perhaps  even  ocean  sources.  The  models  to  be  developed  in  this 
report  deal  strictly  with  surface  water  supplies  in  the  form  of 
rivers  and  reservoirs. 

Many  regions  of  the  world  and  of  the  United  States  in 
particular  depend  upon  rivers  as  their  primary  source  of  fresh 
water  supplies.  Areas  that  depend  upon  this  form  of  water  supply 
are  highly  dependent  upon  rainfall  as  their  source.  Consequently, 
during  periods  of  low  rainfall  or  drought  conditions,  river  flow 
may  be  dangerously  low,  effecting  both  the  amount  and  quality  of 
water  supplied.  Conversely,  during  periods  of  heavy  rainfall,  the 
surrounding  communities  are  dependent  upon  the  river's  ability  to 
remove  the  excess  water  and  to  prevent  serious  and  costly  flooding 
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conditions.  ?or  these  reasons  and  others,  dams  have  been 
constructed  along  existing  rivers  which  back  up  the  water  above  the 
dam,  creating  man-made  reservoirs.  Figure  3-3  shows  a  hypothetical 
river  system  with  two  reservoirs.  The  watershed  for  a  given 
reservoir  is  that  geographical  area  whose  runoff  ultimately  drains 
into  the  reservoir.  In  the  case  shown  the  watershed  for  reservoir 
I  includes  all  of  the  area  upstream  from  the  dam,  whereas  the 
watershed  for  reservoir  2  includes  the  area  between  the  two  dams. 
Reservoir  2  also  receives  water  from  the  releases  of  reservoir  1  . 

Many  users  draw  their  water  directly  from  the  reservoirs 
rather  than  from  the  river.  Since  these  reservoirs  act  as  large 
holding  tanks,  the  3upnly  of  water  can  be  regulated  partially 
through  the  operation  of  the  dams.  The  regulation  of  water  supply 
tends  to  reduce  the  possibility  of  low  water  supply  conditions  by 
storing  up  water  for  approaching  dry  seasons  and  can  act  as  a  flood 
control  system  by  lowering  the  reservoir  level  prior  to  pending 
rainy  seasons. 

Recognizing  that  decisions  mu3t  be  made  periodically 
(daily,  weekly,  monthly,  etc.)  as  to  the  operation  of  a  system  of 
reservoirs,  it  is  logical  to  consider  this  as  a  multiperiod 
decision  problem.  Within  each  period,  decisions  must  be  made  as  to 
how  much  water  to  supply  to  each  of  the  users,  how  much  to  release 
downstream  and  thus  how  much  to  hold  in  the  reservoir  for  use  in 
the  next  period.  The  model  which  will  be  developd  in  thi3  report 


will  provide  this  information  for  the  multiperiod  problem. 


Figure  3-3 

Hypothetical  River  with  Two  Reservoirs 
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3.2.2  The  Single  Period  Model 

In  thi3  section  the  network  model  for  a  single  period  of 
tine  is  constructed.  In  a  later  section  it  is  expanded  to  multiple 
periods.  The  single  period  model  compresses  all  flows  on  a 
particular  facility  (is.,  inflow,  river,  demand,  or  reservoir)  into 
a  single  number  which  represents  the  total  flow  for  the  period. 
Thu3  all  detail  on  flow  variations  within  the  period  are  lost. 
This  discretizing  of  time  is  a  necessary  approximation  to  make  the 
model  of  this  report  computationally  practical.  The  selection  of 
the  time  period  is  an  important  step  in  the  modeling  process. 
Different  applications  might  lead  to  different  selections.  Thus,  a 
time  period  of  a  day  or  even  several  hours  might  be  necessary  for 
the  control  of  a  flood  condition,  while  a  planning  model  for  wafer 
supply  could  use  a  model  with  monthly  or  seasonal  periods. 

First  the  demands  for  water  are  modeled.  For  convenience 
the  models  of  this  report  3how  all  users  drawing  water  directly 
from  the  reservoirs.  All  users  at  a  particular  reservoir  are 
combined  into  a  single  equivalent  user.  It  is  easy  to  enrich  the 
model  for  more  complex  arrangements  by  adding  more  nodes  and  arcs 
along  the  river  reaches.  Demand  is  not  a  fixed  withdrawal  of 
water,  but  rather  i3  measured  for  each  reservoir  by  a  monetary 
benefit  function  for  water  used.  Realistically  one  would  expect 
that  a  benefit  function  would  be  a  concave  function  exhibiting 
decreasing  marginal  return  as  illustrated  in  Figure  3-4.  For  this 
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Figure  3-4 

Total  Benefit  for  Water  Provided  at  Reservoir  1 


Marginal  Benefit  of  Water  Provided  at  Reservoir  1 
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report  a  piecewise  linear  approximation  for  benefit  13  used  with 
the  marginal  return  decreasing  in  steps  as  shown  in  Figure  3-5. 
Although  it  is  recognized  that  providing  a  benefit  function  of  the 
type  described  is  not  an  easy  task,  it  is  required  that  one  be 
estimated  for  each  equivalent  user. 

The  network  model  representing  users  for  each  reservoir  i3 
illustrated  in  Figue  3-6.  Sach  reservoir  i3  represented  by  a  node 
in  the  single  period  model.  For  the  example,  nodes  1  and  2 
represent  reservoirs  1  and  2  respectively.  3ach  user  is  al30 
represented  by  a  node.  A  node  is  also  provided  for  the  ocean.  The 
three  arcs  connecting  each  reservoir  node  with  the  associated  user 
node  represent  the  piecewise  linear  approximation  for  the  benefit 
of  water  to  the  users.  Since  the  network  model  uses  only  cost,  the 
benefits  are  3hown  as  negative  cost3.  Sach  of  these  arcs  has  a 
known  capacity  which  represents  the  step  in  the  piecewise  function. 
Figure  3-5  shows  the  marginal  benefit  for  water  provided  at 
reservoir  1.  Thus,  if  0  to  5  units  are  provided  to  user  I,  the 
benefit  is  330  per  unit.  The  marginal  benefit  for  5  to  1 0  units  is 
320  and  the  marginal  benefit  for  10  to  15  units  is  310.  To 
represent  the  benefit  function  as  a  cost,  one  must  take  the 
negative  of  the  benefit  function.  The  cost  function  thus  formed  is 
convex.  A  negative  slack  external  flow  must  be  provided  at  each 
user  node  to  allow  flows  to  leave  the  network. 

Another  important  aspect  of  the  single  period  model  has  to 
do  with  river  flow.  Rivers  are  represented  by  arcs  between 
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reservoir  nodes  or  from  a  reservoir  node  to  the  ocean  node. 
Conditions  may  exist  whereby  an  especially  low  or  high  flow  would 
be  undesirable.  Consequently,  additional  arcs  can  be  added  to  this 
network  which  would  place  a  premium  on  meeting  certain  low  flow 
conditions  and  a  cost  on  high  flow  conditions,  network  models  allow 
a  lower  bound  on  an  arc  which  forces  the  system  to  supply  at  least 
a  specified  minimum  amount  of  flow  to  that  arc.  However, 
conditions  may  exist  in  which  there  is  not  enough  water  available 
to  provide  even  this  minimum  amount  of  flow,  thus  resulting  in  an 
infeasible  solution.  Another  way  to  handle  thi3  which  is  more 
general  in  nature  and  circumvents  thi3  drawback  is  3hown  in  Figure 
l--7. 

Here,  the  arcs  between  reservoirs  represent  a  piecewise 
linear  function  where  each  arc  has  a  given  capacity.  The  capacity 
of  the  arc  with  a  cost  of  -10  may  represent  the  desired  minimum 
flow  in  the  river  during  the  period.  The  negative  cost  will  tend 
to  provide  that  flow  if  enough  water  is  available  in  the  system  and 
if  other  needs  (also  measured  by  negative  costs)  are  not  more 
important.  The  capacity  of  the  arc  with  zero  cost  would  represent 
the  safe  flow  levels  of  the  river,  between  low  flow  and  flood 
conditions.  The  arc  with  the  high  positive  cost  would  represent 
flood  conditions  and  its  capacity  should  be  set  very  large  (again 
30  as  not  to  create  an  infeasible  solution  in  the  event  of 
extremely  large  quantities  of  water  available).  The  cost  is  the 
penalty  of  allowing  a  flooding  condition.  All  of  these  arcs  would 


have  zero  as  the  lower  bound  on  flow.  Naturally,  additional  arcs 
could  be  used  which  might  be  necessary  to  properly  reflect  the 
results  of  rapidly  increasing  hazardous  and  costly  flood  conditions 
which  would  obviously  not  be  linear  from  the  onset  of  a  flood  to  a 
massive  flood  condition.  These  river  level  arcs  measure  costs  to 
the  current  period.  The  total  cost  curve  of  the  three  arcs 
originating  at  node  1  and  terminating  at  node  2  appears  as  figure 
3-3.  Again,  this  C03t  function  is  convex. 

For  the  single  period  case,  the  water  level  in  the 
reservoirs  may  al30  be  important  to  the  area  for  such  things  as 
wildlife,  sports  and  environmental  issues.  Of  primary  concern 
might  be  the  quality  of  water  if  the  reservoir  is  allowed  to  go  too 
low  and  the  safety  of  local  areas  if  the  reservoir  is  allowed  to 
rise  too  high.  These  concerns  can  be  reflected  in  the  network  by 
using  parallel  arcs  to  represent  each  reservoir.  That  is,  provide 
arcs  to  indicate  minimum,  acceptable,  and  maximum  reservoir  levels 
just  as  was  done  for  the  river  reaches.  This  could  be  done  as 
shown  in  Figure  3-9.  Here,  the  arcs  from  node  1  to  node  la  have 
been  added  to  represent  (-'0)  the  low  condition,  (0)  the  safe  range 
and  (10)  the  high  level  case.  Capacities  on  these  arcs  will 
indicate  the  ranges  over  which  the  given  costs  are  applicable. 
These  arcs  perform  the  very  3ame  function  as  the  river  reach  arcs 
and  their  total  cost  curve  would  be  the  same  form  as  Figure  3-3. 

The  flows  in  the  reservoir  arcs  of  Figure  3-9  represent  the 
water  stored  at  the  end  of  the  period  for  use  in  following  periods. 
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Figure  3-8 

tal  Cost  of  Flow  in  the  River 
from  Reservoir  1  to  Reservoir  2 


>fh  y  a*  ****>»'*■  • 


Iain  factors  couli  be  used  on  the  reservoir  or  river  arcs  to 
represent  losses  due  to  evaporation  or  seepage  (the  gains  would  be 
less  than  one).  For  the  single  period  problem,  water  allocated  to 
the  reservoir  arc3  will  leave  the  network.  Thus  nodes  la  and  2a 
would  have  negative  slack  external  flows.  Since  the  flow  in 
reservoir  arcs  is  limited  by  arc  capacities,  the  slack  external 
flow  should  be  at  least  as  large  as  the  3um  of  the  arcs  entering 
the  node  and  have  a  slack  cost  of  zero.  For  the  multiperiod  case 
nodes  la  and  2a  will  be  nodes  in  the  network  model  of  the  following 
pe riod. 

Ml  that  remains  for  this  single  period  model  is  to  provide 
a  source  of  water  to  the  network.  Inflows  may  be  of  several  types, 
including  runoff  and  ground  seepage  due  to  rainfall,  returns  from 
urban  and  industrial  users  and  imported  water  as  well  as  the 
reservoir  waters  3aved  from  the  previous  period.  All  of  these 
will  be  represented  by  positive  fixed  and  3lack  external  flows  at 
the  nodes  of  the  network.  Fixed  external  flows  are  used  for  inputs 
that  are  not  optional  and  must  be  forced  on  the  network  such  as 
deterministic  runoff,  return  waters  and  reservoir  contents  at  the 
beginning  of  the  period.  Slack  external  flows  can  be  used  for 
optional  inputs  such  as  imported  water. 

Note  that  allowing  inflows  at  nodes  discretizes  the 
locations  of  the  inflows.  Thus  although  runoff  and  irrigation 
return  flows  are  nonpoint  inflows,  they  are  approximated  as  point 
inflows.  The  effects  of  this  approximation  are  diminished  if  more 
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rivar  reaches  (hence  more  nodes)  are  defined. 

The  inflows  to  the  system  caused  by  rainfall  are,  of 
course,  not  known  with  certainty.  It  i3  this  aspect  of  the  system 
model  that  will  receive  the  most  attention  in  the  chapters  to 
follow.  Various  assumptions  will  be  made  about  the  knowledge  of 
water  inputs  and  the  decision  options  available  to  the  controller 
of  the  system.  In  this  chapter,  it  is  assumed  that  all  external 
flows  are  deterministic  and  thus  lenown  with  certainty.  The  values 
chosen  can  be  the  expected  value  derived  from  statistical  analysis 
of  historical  runoff  records  or  they  could  be  specific  historical 
sequences  imposed  on  the  system  to  measure  the  effectiveness  of  the 
system. 

The  entire  two  reservoir  single  period  model  as  a  network 
is  shown  in  Figure  3-10.  All  parameters  shown  previously  were  for 
illustrative  purposes  and  are  not  intended  to  represent  realistic 
values.  Future  network  representations  will  be  of  this  form, 
however,  in  most  cases  the  multiple  arcs  between  nodes  will  be 
3hown  as  a  single  arc  for  clarity. 


3.2.3  The  Multiperiod  Model 

For  the  single  period  system  discussed  above,  it  is  assumed 
that  the  decision  maker  or  system  controller  has  access  to  the 
required  data  to  provide  the  appropriate  measures  of  benefit  or 
cost  for  the  various  network  components  (demand,  river  levels, 
reservoir  levels)  as  well  as  the  rainfall  data.  It  has  also  been 
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assumed  that  these  benefits  can  be  represented  as  piecewise  linear 
functions.  These  benefits  are  strictly  in  terms  of  current  period 
benefits,  that  is,  no  benefits  have  been  specified  beyond  the 
current  period.  For  the  single  period  problem,  the  decision  maker 
observes  the  initial  reservoir  contents,  calculates  an  expected 
inflow  and  solves  the  single  period  model  to  determine  his  optimum 
decisions.  This  process  is  repeated  in  exactly  the  same  way  for 
each  successive  period.  There  is  a  major  drawback  to  this  single 
period  decision  approach.  The  decisions  made  in  one  period  have  a 
direct  effect  on  the  possible  options  available  in  the  next  and, 
in  fact,  in  several  of  the  succeeding  periods.  The  decisions  in  a 
given  period  should  be  made  to  maximize  not  only  the  current 
benefits,  but  also  the  future  benefits.  This  will  be  done  by 
introducing  a  multiperiod  model  that  explicitly  takes  account  of 
the  tradeoffs  between  current  and  future  uses. 

The  multiperiod  model  is  constructed  by  providing  a  single 
period  model  for  each  period  under  consideration.  The  single 
period  models  are  linked  by  the  reservoir  arcs  as  shown  in  the  four 
period  example  of  Figure  3-11.  The  reservoir  arcs  allow  water 
stored  at  the  end  of  one  period  to  be  used  in  the  following 
periods.  The  amounts  of  water  stored,  represented  by  arc  flows, 
are  variables  of  the  optimization.  For  simplicity,  Figure  3-1 1 
only  3hows  one  arc  in  each  of  the  parallel  arc  sets  of  Figure  3-10. 

The  node-arc  structure  of  the  single  period  models  in  the 
combination  are  the  same,  but  the  external  flows,  arc  costs  and  arc 


Figure  3-11 


Four  Period  Model 
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capacities,  which  represent  the  supply  and  demand  of  water,  will 

undoubtably  vary  over  time.  For  instance,  the  example  of  Figure 

3-11  might  represent  a  one  year  time  interval  with  the  four  single 

period  models  representing  the  four  seasons.  Figure  3-12  provides 

a  more  general  schematic  of  a  larger  multiperiod  model.  Here,  the 

single  period  models  are  represented  by  boxes  with  flows  in  period 

t  given  by  the  vector  F,..  Inputs  and  outputs  are  shown  by  the 

vector  quantities  I,,  and  0f.  The  reservoir  arcs  interconnect  the 

periods.  St  i3  a  vector  of  water  quantities  stored  at  the  end  of 

oeriod  t.  Thus,  5  is  the  initial  reservoir  contents  and  3_  is  the 
o  T 

reservoir  contents  at  the  time  horizon.  For  the  deterministic 
model,  all  inflows  must  be  given  (perhaps  in  the  form  of  a 
historical  sequence  of  runoff  data).  Optimization  of  the  network 
model  will  provide  a  policy  for  operation  of  the  system  in  each 
period  3tated  in  terms  of  the  arc  flows. 

One  possible  way  of  utilizing  such  a  solution  in  a 
3tocha3tic  situation,  where  inflows  are  actually  unknown,  i3  to  U3e 
the  solution  as  a  guide  for  policy  in  the  first  period.  Then  when 
the  actual  inflows  are  known  for  the  first  period  along  with  the 
final  reservoir  contents,  the  model  can  be  solved  again  to  obtain  a 
better  solution  for  the  second  period.  The  process  continues  as 
time  progresses  by  solving  the  problem  for  each  new'  period  as  the 
lata  for  the  previous  period  becomes  known. 


3.3  Solution  Procedure  for  the  Pure  Deterministic  Problem 


The  network  models  which  are  used  to  describe  the 
deterministic  water  resource  system  have  the  general  mathematical 
form  of  Model  II  which  is  repeated  here  for  convenience. 

Model  II 

(la) 


Objective : 


.  1  b) 


Minimize  Z 


=±V’ 


Constraints:  Conservation  of  flow  at  each  node: 


f,-r>  ak'k=bi  for  ilin 


k«  M^  ’<«  M. 


(1c) 

Arc  capacity: 


Oif^c^  for  k£  M 

Since  this  problem  is  a  linear  programming  problem,  the 
well  known  methods  of  linear  programming  should,  and  do  provide  a 
solution  procedure.  Thi3  section  describes  in  general  the  primal 
simplex  procedure  as  specialized  to  the  network  flow  problem. 
Computer  solutions  of  network  optimization  problems  are  described 
in  a  book  Ne two rk  Plow  Programming  by  Jensen  and  Barnes  (1990). 
The  procedures  described  in  the  remainder  of  this  report  rest 
heavily  on  the  contents  of  thi3  book.  This  section  is  provided  to 
survey  the  conceptual  ideas  of  the  simplex  technique  applied  to  the 
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generalised  network  minimum  cost  flow  problem. 

5.5-1  The  Primal  Solution 

Svery  linear  program  has  an  optimal  solution  which  is  a 
basic  solution.  For  the  network  problem  a  basic  solution  i3  a 
selection  of  n-1  arcs  (variables)  which  form  an  independent  3et.  A 
selection  forms  an  independent  3et  if  the  columns  from  the 
conservation  of  flow  equations  (lb)  associated  with  the  set  has  a 
nonzero  determinant.  The  basis  for  the  generalized  problem  will 
always  be  a  collection  of  a  single  tree  rooted  at  the  slack  node 
and  zero  or  more  3emi-trees  which  include  a  cycle.  An  example 
problem  is  shown  in  Figure  5-(5.  A  basis  for  this  example  is 
illustrated  in  Figure  5-14.  A  tree  is  a  collection  of  arcs  on 
which  no  cycle  can  be  formed  (neglecting  arc  directions).  A 
semi-tree  will  have  a  single  cycle  with  perhaps  trees  rooted  at 
nodes  on  the  cycle.  Trees  and  3emi-tree3  will  always  be 

represented  by  directing  the  arc3  in  such  a  way  that  there  i3  a 

directed  path  from  the  root  to  every  node.  This  may  necessitate 
reversing  the  direction  of  certain  basic  arcs.  This  is  done  by 
including  mirror  arcs  in  the  tree.  A  mirror  arc  is  given  the  index 
-k  (corresponding  to  the  forward  arc  k) .  While  arc  k  originates 
and  terminates  at  nodes  i  and  .j  respectively,  arc  -k  originates  and 
terminates  at  nodes  .j  and  i  respectively. 

Once  a  basis  is  chosen  the  remaining  arcs  are  called 

nonbasic  arcs.  A  basic  solution  is  formed  by  first  setting  the 


Figure  3-13 

Example  Generalized  Network  Problem 
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flow  on  each  nonbasic  arc  to  either  zero  or  to  the  capacity  of  the 
arc.  The  flows  on  the  basic  arc3  are  then  set  to  values  which 
assure  conservation  of  flow  at  each  node.  Given  flows  for  the 
nonbasic  arcs,  the  flows  for  the  basic  arcs  are  uniquely  determined 
by  the  external  flows  at  the  nodes.  The  basic  solution  may  or  may 
not  be  feasible.  It  is  called  a  feasible  solution  if  the  flow  on 
each  basic  arc  satisfies  the  bounds  on  the  arc: 

0< f,<c,_  for  !c£M„ 

Here  is  the  3et  of  arcs  in  the  basis. 

There  are  many  basic  solutions.  For  each  selection  of  n-1 
independent  arcs  to  form  a  basi3,  there  are  possible  ways  to 
assign  flows  to  the  nonbasic  arcs.  There  may  be  as  many  as  (  ) 
ways  to  choose  the  basic  arcs.  Thus,  an  upper  bound  on  the  number 
of  basic  solutions  is: 

r  m  ■\-3m-n*l 
n-t  )Z 

Linear  programming  theory  tells  us  that  if  a  feasible  solution 
exists,  at  least  one  of  this  large  but  finite  set  will  be  an 
optimum  solution.  It  is  up  to  the  optimization  algorithm  to  find 
and  identify  which  one. 

3-3.2  The  Dual  Solution 

Associated  with  every  linear  programming  problem  is  another 
linear  programming  problem  called  the  dual  problem.  The  dual  of  the 
network  problem  will  not  be  described  here  but  the  dual  variables 
will  be  used  to  check  a  basic  solution  for  optimality  and  to  direct 


60 


the  3earch  for  an  optimum  basic  solution.  The  dual  variables  are 
associated  with  the  nodes  and  are  called  the  node  potentials.  The 
node  potential  for  node  i  is  symbolized  as  7f^. 

For  a  given  basic  network  there  is  a  corresponding  dual 
solution  which  can  be  found  by  requiring  for  each  basic  arc  k(i,.j) 
that  the  following  equality  hold: 

(2) 


/a1(C  for  k(i,.j)€Mg 

If  a  mirror  arc  -k(j,i)  is  in  the  basis  it  i3  required  that: 

(3) 


T^('#>h_k/a_!c)  for  -k  £  Mg 

Assigning  the  parameters  to  the  mirror  arc  in  relation  to  the 
forward  arc  k(i,j)  is  done  as  follows: 

(4) 


h-k=’hk/ak 
’-k"1  /ak 


Squation  (3)  combined  with  equations  (4-)  yields: 

7ri~(7rr\/\)\ 

or 

Tff«r L*\]\ 

which  is  equivalent  to  equation  (2). 

Mote  that  equation  (2)  defines  a  set  of  n-1  linear 
equations  in  n  variables.  Arbitrarily  assigning  zero  as  the 


potential  of  the  slack  node,  the  solution  of  the  equations  then 
yields  the  values  of  the  dual  node  potentials. 
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Now  given  a  basis,  the  primal  solution  (flows)  and  the  dual 
solution  (potentials)  can  be  calculated.  Linear  programming  theory 
provides  the  relations  between  primal  and  dual  solutions  that  can 
be  checked  to  ascertain  optimality.  These  are  as  follows: 

1.  For  each  basic  arc  (k£Mg)  we  have: 

(o)  Primal  feasibility: 

2.  For  each  nonbasic  arc  (k^  M^)  we  have: 

Complementary  slackness: 

(6)  a.  (TT-t-h^/a^  iaplisa  f^'c^ 

(7)  b.  (7f.  *h.  ) / a.  >  7f.  implies  f  *0 

1  :<  .<  J  (C 

(3)  c.  (7f.  *h,  )  /a.»  7 f.  implies  f,  *0  or  2, 

1  fv  !v  ,J  X  !v 

If  given  basic  primal  and  dual  solutions  satisfy  both 
primal  feasibility  and  complementary  slackness  both  solutions  are 
optimal  for  their  respective  problems,  thus  a  test  for  optimality. 

Figure  "5-1 5a  shows  a  flow  solution  for  the  example  problem. 
Figure  3-i5b  shows  the  associated  basic  network  with  node 
potentials  that  satisfy  equation  (3).  It  is  apparent  that  the  flow 
solution  is  basic  and  feasible.  It  only  remains  to  check  the 
complementary  slackness  conditions.  Checking  these  for  each  arc 
reveals  that  the  conditions  are  satisfied  indicating  that  the  flow 
solution  is  optimal. 

3.3.3  Primal  Simplex  Algorithm 

Now  that  there  is  a  procedure  for  checking  optimality,  a 
procedure  is  needed  for  directing  and  carrying  out  the  search  for 


an  optimum  solution.  This  section  describes  a  primal  solution 
procedure  which  is  designed  to  start  from  a  primal  feasible  basic 
solution  and  move  in  a  finite  number  of  iterations  to  an ‘optimum 
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basic  solution. 

First  the  means  to  start  the  procedure  when  no  initial 
basic  feasible  solution  is  known  must  be  provided.  In  general, 
there  is  no  guarantee  that  a  particular  problem  has  a  feasible 
solution.  As  frequently  done  in  general  linear  programming  an 
artifieal  basis  is  used  for  the  network  problem.  Here,  for  each 
node,  an  additional  arc  is  provided  that  connects  the  node  to  the 
slack  node.  The  added  ares  are  called  artifieal  arcs.  These  arcs 
are  generated  according  to  the  following  rules: 

1.  If  bi>0,  create  an  arc  from  i  to  n  with  capacity  b^ 

2.  If  b^O,  create  an  arc  from  n  to  i  with  capacity  -b^ 

3.  For  each  artifieal  arc: 

a.  Assign  the  arc  cost  of  R  where  R  is  a  large 

positive  number 

b.  Assign  a  gain  of  unity 

c.  Assign  a  flow  equal  to  capacity 

4.  Augment  the  arc  set  by  the  artifieal  arcs. 

For  the  artifieal  solution,  all  flows  in  the  original  network  are 
zero  and  all  external  flows  are  carried  on  artifieal  arcs  to  or 
from  the  slack  node.  Conservation  of  flow  is  satisfied  for  all 
nodes,  and  arc  flows  satisfy  flow  bounds.  Thu3,  the  given  solution 
is  feasible  for  the  network  augmented  by  the  artifieal  arcs.  The 
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solution  is  very  expensive,  however,  since  all  flows  pass  through 
the  artifical  arcs  with  marginal  cost  R.  The  algorithm  will 
attempt  to  drive  all  flows  off  of  the  artifical  arcs  and  on  to  the 
arcs  of  the  original  network  to  reduce  the  cost.  If  the  optimal 
solution  has  nonzero  flows  on  any  of  the  artifical  arcs  it  is  clear 
that  there  must  not  be  a  feasible  solution  to  the  original  network 
problem. 

With  an  initial  basic  solution  defined  for  primal  and  dual 
problems,  an  algorithm  that  can  check  for  optimality  i3  required 
and  in  the  event  of  a  nonoptimum  solution,  suggest  a  change  that 
will  bring  the  solution  closer  to  optimality.  This  algorithm  is 
the  primal  basic  simplex  algorithm  which  is  comprised  of  the  3teps 
which  follow: 

1.  Check  each  nonbasic  arc  for  complementary  slackness, 

ie.  , 


if  / Sjc<  tf'j  then  f^c^ 

if  '^i*hk)/ak>  ^  then  f^O 


If 

each 

nonbasic  arc  does  not 

violate  either  of  these 

conditions. 

stop 

,  the 

solution  is  optimal. 

Otherwise,  choose 

an 

arc  to  enter 

the 

basis 

that  violates  one  of 

these  conditions. 

Let 

thi3  be  arc 

*5- 

2.  For  each  arc  in  the  basi3,  find  the  amount  of  flow 
change  in  the  arc  per  unit  of  flow  change  in  arc  k„.  Use  thi3 

It 

information  to  find  the  maximum  flow  change  in  arc  kp  that  will 
cause  the  flow  in  one  of  the  basic  arcs  to  go  to  a  bound  or  cause 
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the  flow  in  arc  k_  to  30  to  it3  opposite  bound.  Choose  the  arc  to 
leave  the  basis  k^  as  the  arc  which  limits  the  flow  change. 

3.  Change  the  flow  in  arc  and  the  basis  arcs  by  the 
amount  found  in  step  2.  If  k^=kP,  return  to  step  1.  Otherwise, 
change  the  basis  tree  by  deleting  arc  k^  and  inserting  arc  kg. 
This  may  require  some  redirection  of  arcs  to  obtain  directed  trees 
and  3emi-trees.  Change  the  node  potentials  to  be  consistent  with 
the  new  basis  network.  Return  to  step  1. 

The  details  of  the  implementation  of  this  algorithm  are 
fairly  complex.  Jensen  and  Barnes  (1980)  provide  complete  details. 
Because  of  the  special  structure  of  the  network  problem  this 
specialized  version  of  the  simplex  algorithm  is  much  more  efficient 
than  more  general  linear  programming  algorithms  applied  to  the 
network  problem. 

* 

3.4  Application  to  the  Guadalupe  River  3asin 

Specific  application  of  this  deterministic  model  was  made 
to  the  Guadalupe  River  Basin  in  Texas.  The  geographical  layout  of 
the  basin  is  as  shown  in  Figure  3-' 6. 

A  proposed  plan  for  the  basin  is  to  expand  the  existing 
reservoir  system  to  include  three  new  reservoirs;  Cloptin  Crossing, 
Cuero  I,  and  Cuero  II.  These  new  reservoirs  along  with  the 
existing  reservoir,  Canyon,  are  felt  necessary  for  meeting  water 
supply  demands  for  the  future.  Future  demands  are  those  projected 
for  the  year  2020  and  are  primarily  for  the  San  Antonio  area,  most 
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of  which  will  be  drawn  off  through  Victoria. 

Monthly  historical  rainfall  and  runoff  data  were  made 
available  by  the  Texas  Department  of  Water  Sesources  for  the  years 
1 925-1 970.  These  data  were  adjusted  to  reflect  runoff  amounts  into 
these  four  reservoirs,  given  that  they  had  existed  during  this  46 
year  period. 

This  four  reservoir  system  i3  shown  in  network  schematic 
form  in  Figure  3-17.  In  contrast  to  the  earlier  networks,  only 
single  arcs  are  3hown  between  the  various  nodes.  This  is  for  the 
3ake  of  clarity  as  multiple  arcs  were  •  sed  for  thi3  example.  This 
network  differs  3lightly  from  the  previous  aodels  in  that  demands 
are  allowed  at  junction  points,  Seguin  and  Victoria,  as  well  as 
from  some  of  the  reservoirs.  This  was  done  to  provide  a  more 
realistic  picture  of  the  true  problem  being  modeled.  AI30,  the 
arcs  between  Cuero  I  and  Cuero  II  running  in  opposite  directions 
imply  an  exchange  capability  due  to  piping  and  pumping  where 
necessary.  This  interchange  was  originally  planned  to  be  an 
equalization  channel  thus  implying  both  reservoirs  would  always  be 
at  the  same  level.  If  modelled  thi3  way,  these  two  reservoirs 
could  be  treated  as  one. 

Twelve  copies  of  the  single  period  model  of  Figure  3-1 7 
were  interconnected  to  form  a  multiperiod  model  with  each  period 
equivalent  to  a  month.  The  multiperiod  model  thus  represents  one 
year  of  operation.  The  operation  of  a  reservoir  system  over  a 
period  of  years  naturally  implies  that  water  available  at  the  end 
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of  this  12  month  network  would  be  available  for  the  13th  month  and 
30  on.  This  is  modeled  by  connecting  the  ending  level  reservoir 
nodee  of  period  12  with  the  beginning  reservoir  nodes  of  period  1, 
resulting  in  a  closed  or  looped  system  of  networks. 

By  summing  the  monthly  inflow  data  and  dividing  by  46,  the 
average  monthly  inflows  were  derived.  These  were  used  as  the 
deterministic  inflows  to  the  system.  AI30  specified  were  reservoir 
capacities,  reservoir  levels  and  demands.  All  data  was  converted 
to  1000  acre  feet  equivalents.  ?or  this  exercise,  demands  were 
distributed  evenly  throughout  the  year.  Thi3  most  likely  would  not 
be  the  true  state  of  nature  for  most  demand  points  and  could  easily 
be  changed  to  reflect  a  more  realistic  demand  profile. 

The  deterministic  solution  showed  that  under  these 
conditions,  enough  water  was  available  on  the  average  to  meet  all 
demands  for  all  periods.  Since  there  was  no  penalty  assessed  for 
releasing  water  to  San  Antonio  3ay,  and  no  reward  for  building  up 
the  levels  of  the  reservoirs,  the  reservoirs  tended  to  be  held  at 
their  minimum  level3.  Given  no  penalty  or  reward  for  doing 
otherwise,  thi3  is  what  one  would  expect  with  known  fixed  inflows 
and  demands.  Specific  node,  arc  and  inflow  parameters  for  this 
twelve  period  model  along  with  the  network  flow  results  3re  shown 
in  the  Appendix  as  the  Guadalupe  River  3asin,  Deterministic  Case. 
A  matrix  generator  was  developed  to  take  the  single  Deriod 
information  and  convert  it  to  a  multiperiod  data  set.  At  this 
point,  the  new  data  3et  could  be  modified  to  reflect  changing  arc 


parameter  choices. 

In  the  next  chapter  a  new  model  is  developed  which  takes 
into  account  the  stochastic  nature  of  runoff  and  the  interaction 
between  reservoirs.  Then  in  Chapter  7  this  new  model  will  be 
applied  to  a  hypothetical  three  reservoir  system  and  finally  to  the 
Guadalupe  River  3asis  four  reservoir  system. 


CHAPTER  IV 


4.  Dynamic  Programming  Solution  Approach 

4.1  The  Decision  Process 

This  chapter  will  be  U3ed  to  develop  and  present  the 
dynamic  programming  solution  approach  to  solve  the  multireservoir, 
multiperiod  stochastic  problem. 

Consideration  of  a  multiperiod  model  implies  that  actions 
taken  in  one  period  effect  not  only  the  current  period  but 
following  ones  as  well.  Decisions  in  the  current  period  must  take 
account  of  the  impact  that  these  decisions  will  have  on  the  periods 
to  come.  To  utilize  this  concept,  a  benefit  function  will  be 
derived  for  each  period  in  the  time  horizon  which  measures  the 
value  of  water  to  be  stored  in  the  system  for  future  use.  Once 
these  benefit  functions  are  known,  the  network  can  be  optimized  for 
a  given  starting  position  with  the  objective  of  maximizing  the 
current  and  future  benefits. 

Realistically,  the  benefit  function  for  any  period  t  should 
reflect  the  characteristic  that  as  more  and  more  water  is  made 
available,  there  i3  a  decreasing  marginal  return  or  benefit  to 
society.  This  implies  that  the  benefit  function  must  be  concave. 
The  model  to  be  developed  requires  that  this  be  the  case. 

Additionally,  for  a  multireservoir  system,  one  would  expect 
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some  interaction  to  occur  between  reservoirs,  particularly  those  in 
close  proximity  to  each  other.  Thus,  for  systems  of  reservoirs,  it 
i3  implied  that  the  future  benefit  of  water  stored  in  a  single 
reservoir  should  be  a  function  of  the  amount  of  water  stored  in 
neighboring  reservoirs. 

A  dynamic  programming  approach  is  used  to  derive  an 
expected  benefit  function  for  water  stored  in  each  period  of  the 
time  horizon.  These  benefits  are  represented  as  a  function  of 
reservoir  contents.  They  are  generally  concave,  nonlinear  and  may 
include  cros3-product  terms  that  represent  the  interaction  between 
reservoirs . 

Because  a  functional  approach  is  used  instead  of  a  discrete 
table,  as.  is  commonly  done,  the  dimensionality  problem  usually 
associated  with  dynamic  programming  is  partially  overcome. 

A  network  flow  programming  algorithm  is  used  to  solve 
3ubproblems  generated  by  dynamic  programming.  Since  the  network 
problems  are  partially  composed  of  nonlinear  objective  functions, 
it  has  been  necessary  to  modify  the  network  algorithms  to  handle 
this  case.  Details  of  how  this  is  done  will  be  presented  in 
Chapter  5. 

The  first  part  of  this  chapter  provides  a  review  of 
deterministic  and  stochastic  dynamic  programming.  This  is  followed 
by  a  brief  description  of  the  multireservoir  network  model  to  be 
used  in  the  dynamic  programming  approach  to  the  multireservoir 
multiperiod  problem.  Finally,  the  overall  algorithm  for  the 
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dynamic  programming  derivation  of  the  desired  benefit  functions  is 
presented . 

4.2  Deterministic  Dynamic  Programming 

Dynamic  programming  i3  a  method  of  solving  an  intricate 
problem  by  decomposing  it  into  a  series  of  stages  (of  time,  space, 
etc.)  and  approaching  the  solution  3tepwise.  This  method,  also 
■cnown  as  Recursive  Optimization,  is  based  on  Bellman's  "principle 
of  optimality''  (Bellman  and  Dreyfus  (1  962))  which  states  that  an 
optimal  set  of  sequential  decisions  has  the  property  that  whatever 
the  first  decision  is,  the  remaining  decisions  mu3t  be  optimal  with 
respect  to  the  outcome  of  the  first  decision. 

In  contrast  to  linear  programming,  there  doe3  not  exist  a 
standard  mathematical  formulation  of  "the”  dynamic  programming 
problem.  Rather,  dynamic  programming  is  a  general  type  of  approach 
to  problem  solving,  and  the  particular  equations  used  must  be 
developed  to  fit  each  individual  situation.  Of  importance  to  any 
dynamic  programming  problem  is  the  identification  of  the  stages, 
the  state  variables,  the  transition  equations,  the  decision  set  and 
the  return  function. 

To  illustrate  the  deterministic  dynamic  programming 
approach  consider  a  multiperiod  single  reservoir  problem.  Here, 
the  stages  would  represent  the  several  time  periods  within  the  time 
horizon.  .4  finite  time  horizon  of  T  periods  is  assumed.  Each 


stage  would  then  represent  a  period  of  time  for  which  a  decision  or 
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set  of  decisions  must  be  made. 

The  decisions  to  be  made  at  each  time  period  may  be  how 
much  water  to  give  up  in  the  form  of  supplying  demands  and 
relsases,  or  how  much  water  to  store  in  the  reservoir  for  the  next 
and  future  time  periods. 

To  make  these  decisions,  one  must  know  the  level  of  water 
3tored  in  the  reservoir  at  the  beginning  of  the  time  period.  The 
level  of  water  in  the  reservoir  is  then  the  state  variable.  The 
ending  level  of  water  in  a  reservoir  will  be  equal  to  its  initial 
level  plus  inflows  minus  outflows.  Let  equal  the  value  of  the 
3tate  variable  at  3tage  t,  i.e.  the  reservoir  level  at  the  end  of 
period  t.  The  value  of  the  state  variable  at  stage  t  is  defined  as 
a  function  of  the  value  at  stage  t-1 : 


;t  *  3t-1  +  it  'dt 


where: 


1  is  the  level  at  the  beginning  of  period  t 
Sf  is  the  level  at  the  end  of  period  t 

i^.  is  the  inflow  in  period  t  (assumed  known  for  the 
deterministic  problem) 

df  is  the  decision  on  outflows  for  period  t 
This  is  called  the  transition  equation  and  it  is  defined 

for  all  t  from  1  through  T.  Its  purpose  is  to  uniquely  define  the 

value  of  the  3tate  variable  at  the  input  to  the  next  stage.  Assume 

for  now  that  the  above  defined  terms  can  only  take  on  the  integer 


values  1,2,... ,9,  and  that  any  combination  of  them  that  forces  3  to 
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be  Leas  then  1  will  take  on  the  value  of  1  or  greater  than  9,  the 
value  9-  These  are  referred  to  as  boundary  conditions.  Hence,  for 
each  stage,  the  state  variable  (water  level)  can  take  on  any  of  9 
values. 

Having  defined  the  stages  and  state  variables  consider  now 
the  possible  decisions.  Let  the  decision  to  be  made  at  each  stage 
equal  the  total  outflow.  Remember  this  may  represent  both  water 
supplied  to  users  and  water  released  downstream.  The  decision  can 
range  from  min(3^_(  +it-9)  to  (St,_1+i^.-1)  since  only  available  water 
can  be  released  and  there  must  be  at  least  one  unit  left.  For  each 
level  of  3,  there  are  at  most  9  possible  transitions  that  can  be 
made  depending  on  the  value  of  i  and  on  the  decision  d.  Thus, 
between  any  two  consecutive  stages  there  are  at  most  31  possible 
paths  to  take. 

Hext,  the  return  function  must  be  defined.  This  represents 
the  immediate  cost  of  making  a  given  decision  starting  from  a  given 
state.  The  immediate  cost  of  making  decision  d^  while  in  state 
3^._1  i3  expressed  as  ,dt) .  Based  on  3ellman's  principle  of 
optimality,  the  total  cost  function  is  the  sum  of  the  immediate 
cost  plus  the  cost  associated  with  making  optimal  decisions  from 
the  new  state  to  the  end  of  the  time  horizon.  This  is  expressed 
as : 

7(St_,,dt)*f(3,) 

where  f(%t)  is  the  cost  of  the  optimum  policy  for  periods  t*1 
through  T.  Taking  the  minimum  of  this  over  the  possible  decisions 
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yields  the  desired  recursive  relationship: 

f(St  ^-MinCcCS,.  ,,dt)+f(St))  for  t-1,...T 
dt 

The  state  variable  3^  represents  the  reservoir  level  at  the 
end  of  period  1 ,  at  the  end  of  period  2,  at  the  end  of  period 
3,  and  so  on.  The  state  variable  3^  is  the  reservoir  level  at  the 
end  of  the  time  horizon.  The  quantity  3g  i3  the  initial  reservoir 
level  ( at  time  0) . 

If  a  value  of  f(S^)  for  every  possible  state  ST  i3  assumed, 
the  recursive  equation  for  t=T  can  be  solved.  'Jote  that  solving 
this  equation  for  the  example  requires  a  minimization  over  as  many 
as  nine  decisions  for  each  state  S?_1 .  There  are  nine  different 
values  of  f(31,_1),  one  for  each  possible  value  of  the  state 
variable  3^  . 

With  fvS,^  )  known,  the  recursive  equation  can  be  solved 
for  t*T-1 .  The  procedure  continues  until  the  recursive  equation  is 
solved  for  t*1  .  At  this  point,  the  solution  i3  complete  except  for 
the  recovery  of  the  optimum  solution. 

The  process  of  recovering  the  optimum  is  called  the 
traceback  procedure.  Let  i1_*(3^._1  )  be  the  optimum  decision  found 
for  state  3^._1  by  solving  the  recursive  equation.  Given  an  initial 
value  of  reservoir  contents  Sg,  the  optimum  decision  for  period  1 
i3  d1*(3g).  The  transition  equation  can  be  used  to  find  the 
optimum  value  of  31  : 

o,  »Sg-di  *(.5g)*i1 
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The  value  of  S,  determines  the  optimum  decision  dp*^  )  which  in 
turn  indicates  the  value  of  Sj.  This  traceback  procedure 
ultimately  leads  to  the  complete  optimum  solution  for  the 
deterministic  problem. 


4.3  Stochastic  Dynamic  Programming 

The  above  example  is  illustrative  of  a  deterministic 
problem  in  that  the  3tate  at  the  next  stage  is  determined  by  the 
inflows  during  the  period  which  are  assumed  known.  In  a  realistic 
situation,  the  inflows  are,  of  course,  not  known  with  certainty 
because  they  depend  on  the  variability  of  nature.  With  stochastic 
dynamic  programming,  it  is  not  necessary  to  assume  a  deterministic 
transition.  Rather,  a  probability  distribution  on  the  transition 
is  defined.  Thus,  let  p( S ^  j  d ^ , S ^ _ 1 )  be  the  probability  of 
transition  to  state  3,.  given  that  decision  d^  is  made  starting  in 


3tate  3. 


Because  3ome  transition  must-be  made,  it  is  required 


all  S 


p(st|dt- Vi )’1 


Define  C(3^._1  ,dt,S<.)  to  be  the  cost  of  starting  in  state  , 
making  decision  d^.  and  ending  in  state  S^.  Now  the  stochastic 
dynamic  programming  recursive  equation  can  be  written.  Let  f(S^) 
be  the  expected  value  of  starting  at  state  3^  (at  the  end  of  period 
t),  traversing  to  the  time  horizon,  and  always  making  the  decison 
which  minimizes  the  expected  cost.  Then: 
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p(sJ'1t,St_i  )(C(St_,  ,dt,St)*f(3t)) 


for  t-1 ,2, . . .  ,T 

This  recursive  equation  is  again  solved  backwards  by  first  assuming 
a  value  for  f ( )  and  then  solving  for  f(S,j  1).  This  allows  the 
solution  for  fCS,^)  to  be  obtained.  The  process  continues  until 
the  value  of  f(Sg)  is  determined.  As  this  equation  is  solved  for 
each  discrete  value  of  5  ^ ,  an  optimum  decision  is  found 

dt,*(Si_1).  This  is  the  decision  that  minimizes  the  expected  cost 
from  period  t  to  the  time  horizon  given  that  the  system  is  in  state 
at  time  t-1  . 

Although  the  optimum  decisions  are  known  for  every  state 
value,  the  optimum  set  of  decisions  for  the  entire  time  horizon 
cannot  be  determined.  The  traceback  operation  previously  defined 
for  deterministic  problems  is  not  applicable  for  the  stochastic 
case.  The  traceback  is  not  applicable  3ince  the  transition  is  not 
certain  at  any  stage,  rather,  it  is  governed  by  a  probability 
distribution.  Although  an  optimum  decision  is  determined  for  each 
state,  only  the  first  decision  is  determined  since  only  Sg  i3 
known.  Stochastic  dynamic  programming  models  a  realistic  decision 
process  in  which  decisions  are  only  made  for  the  current  period  on 
the  basis  of  the  current  state  value. 

It  should  be  stressed  that  stochastic  dynamic  programming 
uses  the  criterion  of  expected  value  of  costs.  There  are  other 
criteria  which  might  be  more  appropriate  3uch  as  stochastic 


dominance  (3arnes  et  al.  (1973)).  It  is  not  clear  however  how 
these  could  be  incorporated  into  an  optimization  algorithm  such  as 
dynamic  programming.  All  the  literature  surveyed  by  the  author 
involving  stochastic  dynamic  programming  utilized  the  expected 
value  criterion,  30  this  criterion  will  be  used  in  this  report. 

This  standard  approach  to  dynamic  programming  has  a  major 
drawback:  the  "curse  of  dimensionality".  As  the  number  of  state 
variables  increases,  the  size  of  the  problem  in  terms  of  both 
computer  storage  and  computation  time  becomes  prohibitive.  In  the 
example  above,  there  was  only  one  state  variable  with  nine  possible 
decisions  at  each  stage.  If  another  state  variable  were  added 
(e.g.  the  level  of  a  second  reservoir) ,  there  would  now  be  31 
unique  states  for  each  stage.  The  number  of  possible  states  at 
each  stage  is  equal  to  the  number  of  levels  raised  to  the  (number 
of  reservoirs)  power.  In  a  five  reservoir  problem  there  would  be 
59,04-9  states.  Problems  in  water  resources  frequently  involve 
systems  of  four  or  more  reservoirs. 

In  the  sections  to  follow,  a  technique  is  developed  which 
replaces  the  discrete  vector  f(Sf)  by  a  single  mathematical 
function  fa  benefit  function)  for  each  period,  and  adopts  a  method 
of  sampling  from  the  distribution  of  inflows.  3oth  of  these 
approaches  greatly  contribute  to  a  reduction  in  the  computational 
requirements  of  dynamic  programming  and  thus  allow  somewhat  larger 


4-. 4  Network  Model  for  Multireservoir  Multiperiod  Problem 

?or  illustrative  purposes  and  for  reference  throughout  thi3 
report  consider  a  three  reservoir  network  as  shown  in  Figure  4-1 . 
This  model  i3  very  3imilar  to  the  two  reservoir  model  of  Figure 
3-10  with  one  exception.  Here,  three  new  arcs  (1,2,3)  and  three 
new  nodes  ^4.,9,12)  have  been  added.  These  new  arcs  will  be  used  to 
represent  the  future  value  of  water  to  the  system.  Their  co3t 
functions  will  include  the  nonlinear  portion  of  the  objective 
function  and  in  most  cases  they  will  be  nonseparable.  It  is  the 
combined  cost  function  of  these  three  nonlinear  arcs  that 
represents  the  future  value  of  water  to  the  system  since  the  flow 
in  these  arcs  represents  water  3tored  for  future  use.  The  cost 
functions  for  these  arcs  will  also  include  linear  terms. 

Because  of  the  lynamic  programming  approach  to  be  used,  it 
is  not  necessary  to  connect  the  ending  reservoir  levels  to  the 
beginning  reservoir  levels  of  the  next  period  as  was  done  in  the 
deterministic  multiperiod  case.  Thus,  all  nonlinear  arcs 
representing  final  storage  could  be  terminated  at  a  sinvle  node. 
For  modeling  and  visual  convenience,  three  nodes  will  be  used 
rather  than  a  single  node. 

The  network  model  of  Figure  4-1  i3  all  that  is  required  for 
the  multiperiod  model.  This  greatly  simplifies  the  data  input 
requirements.  3ecause  of  the  dynamic  programming  approach  this 
seemingly  simple  single  period  model  provides  all  the  required 
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information  for  multiperiod  use.  This  will  become  more  obvious  in 
the  next  few  sections. 

Naturally,  as  the  model  represents  different  periods,  any 
or  all  of  the  network  parameters  can  change  to  reflect  changing 
water  availability  and  requirements  over  time.  As  will  be  seen, 
even  if  it  is  not  desired  to  change  the  network  parameters  in  the 
linear  part  of  the  network,  the  cost  parameters  of  the  nonlinear 
arcs  will  necessarily  change  from  period  to  period  due  to  the 
changing  inflow  parameters  causing  flow  changes  in  these  arcs.  A 
change  of  flow  in  any  of  the  nonlinear  arcs  can  cause  a  change  in 
the  cost  assigned  to  other  nonlinear  arcs  if  the  cross  product 
terms  are  nonzero. 
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4.5  The  Dynamic  Programming  Algorithm  for  Deriving  the  Benefit 
Functions  for  the  Multiperiod  Multireservoir  Problem 

Figure  4-2  represents  the  multiperiod  model  for  the 
stochastic  problem.  Each  box  is  a  single  period  model  of  the  type 
shown  in  Figure  4-1.  The  notation  in  this  figure  is  as  follows: 

T  *  number  of  periods  of  the  analysis  or  time  horizon 
R  *  number  of  reservoirs,  r*1,...,R 


( i. . ,i_. , . . . , i  ,).  This  is  a  vector  of  runoff 
i t  it  rt 

inflows  to  the  reservoirs  during  period  t.  Assume  that 


I  is  a  random  vector  from  a  known  distribution.  The 
x 

parameters  of  the  distribution  or  the  distribution 

itself  may  differ  from  period  to  period.  The  inflows 

for  two  different  periods  are  assumed  to  be  independent 

random  variables.  It  is  assumed  that  these  are  the 

only  external  flows  into  the  system  except  the  initial 

reservoir  contents  of  period  1 . 

S.  »  (s, , ,s_. , . . . ,s  .).  This  is  a  vector  of  reservoir 

t  7  x  rx 

contents  at  the  end  of  period  t.  Mote  that  it  also 
describes  the  initial  reservoir  contents  of  period  t*1 . 


\  *  ^°1 t’°2t’ * ' ' ,0rt^ ' 

during  period  t. 


This  is  a  vector  of  outflows 


.).  This  is  a  vector  of  arc  flows  in 
i  x  <sx  mx 

the  single  period  network  for  time  t.  The  vectors  3 


3 
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and  Of  ar9  also  included  as  flows  on  arcs. 

The  inflows  to  the  network  in  period  t  are: 

3t-1+It 

They  are  represented  as  positive  external  flows  at  the  reservoir 
nodes.  The  outflows  of  the  system  are  flows  to  the  demand  nodes 
and  the  end  of  the  period  reservoir  contents  nodes.  For  the  single 
period  model  of  Figure  4-1  it  is  assumed  that  inflows  will  be  known 
before  the  flow  decisions  are  made.  In  reality,  flows  are 
continually  adjusted  throughout  a  period  as  the  inflows  are 
revealed  by  the  passage  of  time.  For  this  discrete  model  the 
length  of  the  time  period  may  be  adjusted  to  allow  any  desired 
degree  of  accuracy  however  the  distribution  of  flows  within  a  time 
period  are  assumed  to  be  instantaneous. 

The  multiperiod  model  is  solved  with  stochastic  dynamic 
programming.  The  flow  chart  of  Figure  4-3  represents  the  overall 
dynamic  programming  process  used  for  deriving  the  benefit  functions 
for  all  t.  The  periods  referred  to  in  this  flow  chart  represent 
the  T  periods  of  Figure  4-2. 

Before  describing  the  algorithm  in  detail  some  additional 
notation  is  necessary. 

Let: 

TO  *  'lumber  of  discretizations  of  reservoir  contents 

TKp  ■  Capacity  of  reservoir  r 

RL  *  Vector  of  reservoir  levels  for  discretization  as  a 


Figure  4-3 

Dynamic  Programming  Approach 


87 


percent  of  CfC?.  There  will  be  NN  components  of  RL. 

ZZ  *  Value  of  the  random  number 

»  Number  of  level  combinations  given  R  reservoirs. 
Thus, 

Lr  -  NN* 

K  *  Number  of  random  observations  per  level  combination 
of  the  vector  If  taken  from  the  distribution  of  1^ 

3P(3t)  =  Expected  benefit  for  water  stored  at  the  end  of 
period  t  for  t  =  1 . T 

.1^)  *  Model  response  for  level  combination  when 

the  random  inflow  is  1^ 

Using  the  above  notation  30me  sample  values  will  be  assumed  for 
purposes  of  illustrating  the  overall  process.  For  this  example, 
let: 

R  *  2  reservoirs 
NN  *  3 

UK,  -  20,  CK2  *  40 
RL  -  (.5.  .7,  .9) 

K  *  10 

Bf(3t)  -  20f ,  *  25f2  -  ff  -  f§  -.45f,f2 

Then: 

Ljj  ■  3^  *  9  level  combinations 

These  9  level  combinations  will  be  represented  by  the  vector  St_,  » 
C •  32t ) "  possibilities  are  listed  in  Table  4.-1. 
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Table  4.-1 

Level  Combinations  for  the  Example 


91t 

32t 

10 

20 

10 

23 

10 

36 

14 

20 

14 

23 

14 

36 

13 

20 

13 

23 

13 

36 

Thi3  example  will  be  continued  later. 

The  steps  of  the  algorithm  corresponding  to  the  boxes  of 
the  flow  chart  in  Figure  4-3  are  as  follows: 

1 .  Read  the  network  and  reservoir  data.  This  includes  reading 
all  of  the  network  arc  parameters  and  all  reservoir  data 
(TO,  CKr,  y,  RL,  T,  etc.).  Also  read  the  runoff  parameters 
for  period  T  for  each  reservoir.  These  define  the 
distribution  of  random  inflows  I  Set  t  *  T. 

2.  Read  SF( 3^ ) :  This  is  the  assumed  functional  representation 
for  the  future  value  of  water  at  the  end  of  period  T.  This 
information  relates  to  the  cost  parameters  to  be  assigned  to 


;  i 

i  1 


* 

•t 


the  nonlinear  arcs  for  period  T.  For  the  two  reservoir 
example,  read  two  linear  coefficients  two  coefficients  for 


K 


89 

the  squared  terms  and  one  coefficient  for  the  cross  product. 

3.  Draw  a  random  observation  from  the  distribution  If . 

4.  Select  a  level  combination.  From  the  combinations, 

select  one  that  has  not  been  evaluated.  This  level 
combination  is  the  vector  since  it  represents  a  given 

initial  set  of  reservoir  contents  for  period  t.  Go  to  5. 

5.  Solve  the  network.  Since  3F(S  ),  and  1^  are  known,  the 

single  period  network  problem  is  deterministic. 

Let  the  optimal  solution  to  the  network  be  given  as: 

Y(St_, ,It)  -  Min(HFt  -  B?(St)) 
subject  to  conservation  of  flow  and  bounding  constraints. 

This  is  not  a  linear  problem  since  BFCS^)  is  nonlinear  and 

3^  contains  variables  of  the  flow  problem.  The  problem  is  a 

nonlinear  network  flow  problem.  The  solution  procedure  for  « 

this  nonlinear  network  will  be  discussed  in  Chapter  5-  To 

further  simplify  the  notation,  let: 

Ti,w  -  » -  1 . K  ; 

This  is  equivalent  to  stating  that  there  are  now  i  different 

i 

level  combinations  and  for  each  of  these  level  combinations  ? 

K  random  draws  will  be  made.  Thus,  Y.  represents  the  J 

value  of  the  network  optimal  solution  for  the  ith  level  and  • 

j 

the  wth  draw.  Go  to  6 

6.  If  all  level  combinations  have  not  been  evaluated,  go  to  4. 

Otherwise,  L-  values  of  Y.  will  have  been  generated  for 

“  i  .w  j 

i 

4 

5 
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thi3  choice  of  v,  one  for  each  level  combination.  For  each 
of  these  level  combinations  the  3ame  random  inflows,  I^, 
will  have  been  added.  (Jo  to  7 . 

7.  If  all  random  draws  ha^e  not  been  made,  go  to  3.  Otherwise, 
for  each  random  draw  all  level  combinations  will  have  been 
evaluated.  A  total  of  L^  times  K  (90  for  this  example) 
optimal  solutions  will  have  been  generated.  Sach  random 
inflow  I  was  applied  to  all  L^  combinations.  This  aspect 
of  the  algorithm  is  discussed  more  fully  in  section  4-. 6.  Go 
to  9. 

9.  Average  the  values  obtained  from  the  observations  for  each 
level  combination.  So  to  9. 

9.  Least  squares  regression.  Perform  a  least  squares 

regression  using  the  averages  of  the  observations  as  the 
dependent  variable.  This  will  be  discussed  more  fully  in 

section  4.7.  Go  to  10. 

10.  Move  back  one  time  period:  Let  t  *  t-1 .  Go  to  11. 

11.  If  t  *  zero,  STOP.  Otherwise,  go  to  11. 

12.  Adjust  network  parameters.  This  requires  adjusting  the  cost 
parameters  for  the  nonlinear  arcs.  These  new  parameters 
will  be  the  coefficents  of  the  derived  benefit  function.  Go 
to  12. 

13.  Read  runoff  parameters  for  period  t-1.  Go  to  3. 

Once  BF(St_1)  is  available,  then  BF(S„_?)  is  computed  in  a 


1 
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like  manner.  The  process  continues  until  BF(51 )  is  evaluated. 
This  approach  is  different  than  the  classical  dynamic  programming 
approach  in  the  following  ways: 

1 .  The  benefit  function  is  represented  as  a  continuous 
mathematical  function  rather  than  for  discrete  values  of  the 
state  variables. 

2.  The  recursive  equation  is  solved  using  a  Monte  Carlo 
sampling  approach  rather  than  using  transition 
probabilities. 

3.  The  optimum  decisions  are  found  using  a  network  flow 
algorithm  rather  than  a  discrete  search  over  a  finite  set  of 
decisions. 

The  solution  procedure  of  the  network  optimization 
algorithm  requires  that  the  form  of  the  benefit  function  be 
specified.  Thi3  is  necessary  since  a  different  solution  technique 
would  be  required  for  solving  the  nonlinear  form  of  the  benefit 
function.  While  any  convex  form  could  be  applied,  the  next  section 
discusses  the  rational  for  using  a  quadratic  to  represent  these 
benefit  functions. 

4.6  Sampling  From  the  Distribution  of  Inflows 

An  important  aspect  of  the  algorithm  is  the  derivation  of 
the  expected  benefit  by  sampling  from  the  distribution  of  inflows. 
The  runoff  distributions  provided  as  3n  input  to  the  procedure  may 
take  many  forms.  These  may  include  normal,  log-normal,  log-normal 
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Pearson  Type  III,  exponential,  historical,  etc.  The  computer 
program  currently  has  three  options  available,  normal,  log-normal 
and  historical.  Other  inflow  distributions  could  easily  be 
incorpo rated . 

As  an  example  of  how  these  inflow  parameters  are  used, 
consider  again  the  example  problem  used  earlier.  If  the 
distribution  of  runoff  is  assumed  to  be  normal,  the  mean  and 
standard  deviation  for  inflows  will  be  read  for  each  reservoir  in 
box  1  of  Figure  4-3.  Let  the  mean  and  standard  deviation  for  the 
two  reservoirs  be  as  follows: 

Reservoir  1  Mean  =  5,  Standard  Deviation  =  2 
Reservoir  2  Mean  =  S,  Standard  Deviation  =  3 
Mext,  in  box  3,  a  random  number  ZZ  will  be  generated  from  a  normal 
distribution  with  a  mean  of  zero  and  a  standard  deviation  of  1 .  ZZ 
is  used  in  conjunction  with  the  mean  and  standard  deviation  for 
each  reservoir.  Thus,  the  total  inflow  for  each  reservoir  will  be 
equal  to: 

It  =  Mean  *  ZZ  *  Std.  Dev. 

If  I^  is  less  than  or  equal  to  zero,  a  zero  inflow  is  assumed.  For 
example,  if  ZZ  3  . 5 ,  i1 1  *  5  *  1  *6,  and  i2t  *  3  +  1  .5  *  9.5. 

This  implies  that  for  the  model  developed,  perfect 
correlation  of  inflows  between  the  reservoirs  in  a  basin  is 
assumed.  This  is  not  a  requirement  of  the  algorithm  and  other 
assumptions  could  be  made. 

In  addition  to  applying  the  random  number  ZZ  to  all 


reservoirs,  this  same  random  inflow  will  be  applied  to  all 
reservoirs  for  all  level  combinations.  For  each  of  the  9  level 
combinations  of  the  example,  add  6  and  9-5  units  to  the  initial 
contents  for  reservoirs  1  and  2  respectively.  This  will  result  in 
total  water  available  for  level  combination  1  of  16  for  reservoir  1 
and  29-5  for  reservoir  2.  For  the  second  level  combination,  these 
values  would  be  16  and  31.5,  etc.  This  process  is  used  to  assure 
convexity  of  the  response  surface.  This  assures  that  if  the  first 
level  combination  is  overestimated,  all  of  the  remaining  level 
combinations  will  also  be  overestimated.  To  minimize  thi3  over  (or 
under)  estimation  it  is  required  that  a  sufficient  number  of  random 
draws  be  made  to  make  this  error  negligible.  This  will  be 
discussed  more  fully  in  Chapter  6. 

4.'7  Least  Squares  Regression 

As  indicated  by  the  algorithm  a  least  squares  regression  is 
performed  on  the  mean  responses  of  the  L^  level  combinatons.  If 
all  level  combinations  have  been  evaluated,  there  will  exist  a 
matrix  of  optimum  solution  values  referred  to  here  as  the  response 
matrix  as  illustrated  in  Figure  4-4-.  The  rows  of  the  response 
matrix  represent  the  level  combinations  and  the  K  columns 
represent  the  network  optimal  solutions,  Y.  for  each  of  the  K 
draws.  Preceding  this  response  matrix  is  another  matrix  which  will 
be  referred  to  as  the  design  matrix.  This  matrix  will  have  rows 
and  9  columns,  where  U  depends  upon  the  number  of  terms  i.n  khe 
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benefit  function  BF(S1._1).  If  it  is  desired  to  fit  a  linear 
benefit  function  to  this  data  in  terms  of  the  R  reservoirs,'  U  *  R  + 
1,  (R  linear  terms  plus  a  constant).  To  represent  the  benefit 
function  as  a  quadratic,  linear,  second  order  and  interactive  terms 

p 

are  required.  Thus,  for  the  quadratic,  U  *  2  *  R  *  +  1 .  Thus, 
the  U  terms  in  each  row  represent  the  desired  form  of  the  benefit 
function.  For  the  example,  U  *  6. 

To  fit  the  desired  form  of  the  benefit  function  to  this 
data,  a  least  squares  regression  i3  performed.  Here,  the  design 
matrix  has  as  its  first  column,  a  vector  of  ones.  This  is 
necessary  to  account  for  the  constant  term.  Accordingly,  this 
design  matrix  represents  the  independent  variables  for  he 
regression  analysis.  As  the  dependent  variable,  the  means  of  the 
rows  of  the  response  matrix  will  be  used.  Thus,  the  dependent 
variable  for  row  i  i3: 


Y. 

i 


for  all  i 


By  fitting  this  data  to  the  selected  design,  BFCS^)  is  derived. 


1.8  The  Quadratic  Benefit  Function 

The  usual  procedure  of  discrete  dynamic  programming  is  to 
store  the  return  function  f(3j.)  in  computer  memory  for  a  large 
number  of  discrete  values  of  the  state  variables  S,..  Vhen  is  a 
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multidimensional  vector  (one  dimension  for  each  reservoir)  the 
number  of  possible  discrete  states  can  be  very  large.  Because 
rapid  access  computer  memory  is  finite  in  3ize,  this  limits  the 
number  of  state  variables  that  can  be  handled  at  each  stage. 
Three  state  variables  is  frequently  described  as  a  practical 
limit. 

In  the  method  presented  here,  this  storage  problem  is 
overcome  by  fitting  a  quadratic  function  to  BF(St)  and  storing  only 
the  coefficients  of  the  quadratic. 

k  quadratic  form  has  been  chosen  for  the  following  reasons: 


1 .  It  can  exhibit  the  concave  shape  expected  for  the 
benefit  function,  (convex  cost  function) 

2.  It  can  represent  nonseparable  interactions  between 
reservoirs  with  cross  product  terms. 

3.  It  is  easy  to  store  in  a  computer. 

4.  It  is  computationally  convenient  in  the  network  models 
which  arise  in  the  solution  procedure. 

If  has  n  dimensions,  the  number  of  terms  in  a  quadratic 
is : 

■*r(n'*3n)  +  1 

For  the  three  reservoir  case  the  full  three  variable  quadratic 
would  have  10  terms  as  shown  here: 


►B5f22*B6f32+B?f 1 f2+Bgf , f3*Bgf?f3 


p 

Storing  these  10  coefficients  is  much  easier  than  storing  the  SB 
discrete  state  vectors. 

Also  beneficial  is  the  fact  that  the  function  is  defined 
for  continuous  St  rather  than  discrete  values.  This  means  that 
once  such  a  function  has  been  derived  (albeit  from  an  approximation 
to  the  discrete  representation)  decisions  can  be  made  by  a  one  time 
solution  of  the  network  using  the  observed  reservoir  levels.  These 
reservoir  levels  need  not  be  equated  or  rounded  to  the  nearest 
discrete  level. 

It  is  clear  that  the  true  expected  benefit  function  is  not 
in  reality  a  quadratic  function.  A  better  fit  to  the  observed  data 
■night  be  obtained  with  a  more  complex  model.  This  research  has 
limited  consideration  to  the  quadratic  because  of  the  reasons  noted 
above.  The  dynamic  programming  methodology  however  is  not  limited 
to  this  case.  Indeed  the  data  could  be  fit  to  any  model  and  the 
recursive  procedure  is  independent  of  the  model.  The  network  flow 
solution  procedure  is  limited  to  the  quadratic  case,  however  it 
probably  would  be  possible  to  derive  a  more  general  procedure  along 
the  lines  of  nonlinear  algorithms  for  pure  network  flow  problems 
'(Luenberger  (1965)).  This  was  outside  the  scope  of  this  research. 

The  next  chapter  will  present  the  methodology  for  solution 
of  the  nonlinear  (quadratic)  network.  The  method  for  obtaining 
data  for  use  in  the  least  squares  program  will  be  fully  explained 
in  Chapter  6. 


Chapter  V 

5.  Solution  of  Nonlnear  Network  Problems 
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5.1  Introduction 

In  the  performance  of  the  dynamic  programming  algorithm  it 
ie  necessary  to  solve  network  problems  with  convex,  quadratic, 
nonseparable  arc  costs.  In  Figure  4-1  the  flow  on  arcs  1 ,  2  and  3 
represent  water  stored  in  reservoirs  for  future  use.  The  benefit 
function  for  this  water  stored  is  represented  as  a  concave, 
quadratic  function.  As  an  example  of  a  quadratic  concave  benefit 
function,  BF(  ft  ,  f2<  f.^)  is  assumed  to  be: 

59f ,  *46f 2>39f3-.86f  1 2-. 53f  . 52f  . 64f  1  fg-. 40f ,  ty. 68f 2f ^ 

where  f1  ,  f2  and  f^  are  the  flows  on  arcs  1 ,  2  and  3  respectively. 
This  benefit  function  will  be  used  later  as  the  period  T  assumed 
benefit  function  for  many  of  the  example  problems  of  Chapter  7. 
Since  the  algorithm  operates  on  costs  the  negative  of  the  benefit 
function  is  used  to  obtain  the  cost  function.  Thus: 

C(ft,f2,f3)  »  -BF(f1tf2,f3) 

This  is  a  convex  cost  function.  The  network  flow  programming 
algorithms  which  have  appeared  in  the  literature  (Ali  at  al. 
(1978),  Cooper  and  Kennington  (1977),  Dembo  and  Klincevicz  (1979), 
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Florian  (1977),  Prank  and  Wolfe  0  956),  Helgason  and  Kennington 
(1973)  and  Klincewicz  (1979)  have  been  designed  for  linear  or 
convex  problems,  but  not  for  generalized  problems  with  nonseparable 
objective  functions.  This  chapter  provides  details  on  the 
theoretical  development  of  an  algorithm  to  handle  quadratic 
nonseparable  objective  functions.  The  procedure  is  designed  for 
generalized  networks  (i.e.  networks  with  gains)  with  convex, 
quadratic,  nonseparable  objective  functions  and  has  been  coded  in 
Fortran  for  the  CDC  computer. 


5.2  Problem  Statement 

Consider  a  network  problem  defined  as  in  Chapter  III  with 
the  added  stipulation  that  a  subset  of  the  arcs,  ,  have  nonlinear 
arc  costs.  This  subset  is  included  in  the  set  of  all  arcs  M.  The 
linear  cost  coefficients  are  described  by  the  vector  H  for  all 
arcs.  A  matrix  will  be  used  to  define  the  nonlinear  component  of 
cost. 


The  cost  of  nonlinear  arcs  is  assumed  to  be  a  quadratic 

function  of  arc  flows.  Let  F„  be  the  vector  of  flows  in  the 

N 

nonlinear  arcs  and  Z„  be  the  nonlinear  cost  contribution  of  the 
nonlinear  arcs.  Then: 


Q  FN 


where  Q  is  a  symmetric  matrix  which  is  positive  definite  and  the  T 
represents  the  transpose  of  a  vector. 
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For  the  example  problem: 


This  is  a  .  positive  definite  matrix,  30  Z  is  a  convex  function. 
This  Q  matrix  is  equivalent  to  the  Hessian  of  CCf^f^.f-j)  divided 
by  2.  This  notation  is  common  in  the  literature  and  Q  is  referred 
to  as  the  quadratic  matrix. 

The  total  cost  for  the  system  flow  is: 

Z-HF+fJ  Q  Fn 

The  nonlinear  arcs  are  also  represented  in  the  flow  vector 
F  so  that  linear  costs  can  also  be  associated  with  the  arcs. 

Define  H'  to  be  the  vector  of  first  derivatives  of  the  arc 
costs.  Of  course,  for  the  linear  arcs: 

Eq.  (1): 

hV\ 

For  nonlinear  arcs: 

Eq.  (2): 

h  k"V29kFN 

where  Q,  is  the  kth  row  of  Q.  In  this  chapter,  Q  will  be 
subscripted  with  k  to  indicate  a  specific  row  in  the  Q  matrix. 
When  Q  is  not  subscripted,  it  will  refer  to  the  entire  matrix.  In 
Chapter  6,  the  matrix  Q  will  be  subscripted  as  Q  when  it  is 
desired  to  identify  it  with  a  specific  time  period. 
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The  algorithm  haa  to  find  a  minimum  coat  aolution  for  the 
flow  vector  P  (which  includea  the  nonlinear  flowa  F„).  The 

IN 

nonlinear  network  model  ia: 


These  are  the  aame  network  conatrainta  aa  for  the  linear  model.  A 


primal  approach  ia  uaed  in  which  an  initial  basic  feasible  aolution 
ia  defined.  This  aolution  describes  a  basis  network.  A  basis 


network  for  the  purs  network  ia  a  set  of  n-1  area  which  fora  a  tree 
rooted  at  the  alack  node  and  having  a  directed  path  from  the  alack 
node  to  all  the  nodes  of  the  network.  For  the  generalized 
network,  the  basis  network  may  consist  of  several  components,  one 
of  which  ia  a  tree  rooted  at  the  alack  node  and  the  other  are  trees 


rooted  at  cycles.  A  component  consists  of  a  set  of  nodes  such  that 
there  is  a  path  between  every  pair  of  nodes  in  the  set. 

The  solution  procedure  for  the  nonlinear  network  parallels 
that  of  the  primal  linear  solution  algorithm  of  Chapter  3.3.  This 
algorithm  is  restated  here. 

1.  Check  each  nonbasic  arc  for  complementary  slackness, 

if  {IT.  *  h,')/a,  <  7f .  then  f,  ■  c, 
l  k  k  j  k  k 


if  ( Tf  ^  *  h^)/a^  >'^j  then  ^  a  ^ 

If  each  nonbasic  arc  does  not  violate  either  of  these  conditions, 
stop,  the  solution  is  optimal.  Otherwise,  choose  an  arc  to  enter 
the  basis  that  violates  one  of  these  conditions.  Let  this  be  arc 
kg*  Note  that  h^  i3  used  in  these  conditions  for  the  nonlinear 
problem. 

2.  For  each  arc  in  the  basis,  find  the  amount  of  flow 
change  in  the  arc  per  unit  of  flow  change  in  arc  k_.  lJse  this 
information  to  find  the  maximum  flow  change  in  arc  k_  that  will 
cause  the  flow  in  one  of  the  basic  arcs  to  go  to  a  bound  or  cause 
the  flow  in  arc  k^  to  go  to  its  opposite  bound.  Choose  the  arc  to 
leave  the  basis,  k^,  as  the  arc  which  limits  the  flow. 

3.  Change  the  flow  in  arc  k^  and  the  basis  arcs  by  the 
amount  found  in  step  2.  If  k^  «  kg,  return  to  step  1.  Otherwise, 
change  the  basis  tree  by  deleting  arc  k^  and  inserting  arc  kg. 
Change  the  node  potentials  to  be  consistent  with  the  new  basis 
network.  Return  to  step  1. 

For  a  linear  problem  the  iterative  step  of  the  primal 
method  checks  all  nonbasic  arcs  for  optimality.  This  is  the  test 
for  complementary  slackness  (Chapter  III)  which  evaluates: 

dk  m7fi  *  K 
■3  3  3  3  JE 


where  the  subscript  3  refers  to  the  entering  arc.  If  d.  is  less 

3 

than  zero  and  f  is  zero,  then  the  network  is  not  optimal  and  arc 
a 

kj,  will  attempt  to  enter  the  basis,  d  is  interpreted  a3  the 
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change  in  the  objective  function  per  unit  change  of  flow  in  arc  k^. 

Since  d  ia  negative,  and  the  model  is  trying  to  minimize  costa, 
kB 

it  would  like  to  put  as  much  flow  on  this  arc  aa  possible.  For 

each  unit  of  flow  added  to  arc  kp,  the  objective  function  will 

change  by  d,  •  In  this  case  the  objective  function  will  be  reduced 

since  d.  <0.  The  other  nonoptimal  condition  is  where  d,  >0  and 
""E  E 

f  »  c,  .  In  this  case  the  algorithm  will  remove  flow  from  arc 
kE  *S 

k_,  and  for  each  unit  removed,  the  objective  function  will  be 

decreased  by  the  amount  d  .  Because  of  the  conservation  of  flow 

*E 

conditions,  changing  flow  on  arc  k_  means  adjusting  the  flows  on 
other  basic  arcs.  The  final  result  is  that  eventually  one  of  these 
arcs,  either  k,  or  a  basic  arc,  will  reach  its  bound  on  flow  and  no 
further  improvement  can  be  made.  If  there  is  enough  slack  in  the 
network  to  allow  arc  k„  to  receive  as  much  flow  as  it  can  handle 

(i.e.  f  goes  to  one  of  its  bounds)  then  arc  k_  i3  not  allowed  to 

Zt 

enter  the  basis  rather  it  becomes  nonbasic  at  the  opposite  bound. 
However,  if  there  is  a  basic  arc  whose  flow  goes  to  one  of  its 

bounds  before  the  flow  on  the  entering  arc  reaches  a  bound  thi3 

basic  arc  will  leave  the  basis  and  arc  k„  will  enter. 

it 

There  are  two  principle  differences  between  the  nonlinear 
and  linear  networks  which  require  special  consideration  in  the 
algorithm.  The  first  is  that  as  flow  changes  on  the  nonlinear 
arcs,  their  marginal  costs,  h,^,  will  also  change.  These  costs  are 
a  linear  function  of  the  flows  and  are  not  constants.  Thus  a  flow 
change  on  the  basic  arcs  and  arc  k„  may  simultaneously  change  the 

Zj 
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dual  variable,  ''  ,  and  the  value  of  h'  .  This  results  in  the 
second  difference.  A  flow  smaller  than  the  one  required  to  cause 
either  k„  or  a  basic  arc  to  reach  a  bound  may  cause  the  entering 

£i 

arc  to  satisfy  complementary  slackness.  Thus,  it  may  be  that  a 

nonlinear  arc  will  remain  nonbasic,  even  though  the  flow  on  the 

nonbasic  arc  is  between  its  bounds.  For  the  linear  problem  each 

nonbasic  arc  will  have  flow  at  zero  or  capacity,  while  for  the 

nonlinear  problem  a  nonbasic  flow  may  be  between  zero  and  capacity. 

These  differences  -will  be  illustrated  below  by  first  demonstrating 

the  algorithm  with  a  linear  network  and  then  showing  the 

differences  with  a  nonlinear  network. 

Assume  that  in  step  1  of  the  primal  simplex  algorithm,  a 

nonbasic  arc  is  discovered  which  violates  complementary  slackness. 

For  this  discussion,  let  d,  <  0.  In  the  linear  network  this  would 

E 

imply  that  f  ■  0.  The  networks  of  Figure  5-1  illustrate  the 
KS 

steps  of  the  primal  simplex  algorithm  for  a  pure  linear  network 
(a^  *  l).  The  basis  for  Figure  5-1a  is  made  up  of  arcs  (2,3,4-). 

The  node  potentials  are  determined  by  the  basic  arcs.  Thus,  7! ^  * 
"^♦h2*0*2»2.  ^  “  2  ♦  3  -  5  and  ^  ^  *  h^  *  5 

*  1  *  6.  Step  1  of  the  algorithm  checks  9ach  nonbasic  arc  for 

complementary  slackness.  The  nonbasic  arcs  are  arcs  1  and  5.  For 
arc  1  ,  d1  -  ^  ♦  h,  -  -  -3  and  d^  -  *  h^  -  Tf  x  -  -1  .  3oth 

arcs  1  and  5  have  d  <  0.  Arc  5,  however,  is  at  capacity  (f_  *  c_ ) 
which  means  that  complementary  slackness  is  satisfied.  Arc  1  has 
zero  flow  and  a  negative  d^  implies  a  nonoptimal  condition.  Arc  1 
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is  identified  as  the  entering  arc,  kg  *  1  . 

Step  2  of  the  algorithm  finds  the  arc  which  will  limit  the 

flow  change  as  flow  is  changed  on  arc  k^.  As  flow  is  added  to  arc 

1 ,  flow  must  be  decreased  on  basic  arcs  2  and  3  to  maintain 
conservation  of  flow.  A  decrease  of  1  unit  on  arc  3  results  in  a 

flow  of  zero  on  this  arc.  This  is  the  limiting  arc  and  k^  *  3. 

Step  3  of  the  algorithm  changes  the  flows  according  to  the  step  2 
results  and  updates  the  node  potentials.  Figure  5-1 b  shows  the 
results  of  this  iteration  of  the  algorithm.  The  basic  arcs  are  now 
(1,2,4).  The  total  cost  for  the  network  flows  went  from  16  to  13 

for  this  iteration.  1,  is  the  change  in  total  cost  for  a  unit 

a 

change  of  flow  in  arc  k_,.  d  was  -3  and  one  unit  of  flow  was 

2,  kg 

added  to  arc  k^  *  1  .  The  total  cost  changed  by  16  -  3  *  13.  The 
algorithm  now  returns  to  step  1  to  check  for  optimality. 

The  basis  for  Figure  5-1 b  includes  arcs  (1,2,4). 
Svaluating  the  nonbasic  arcs  (3,5)  results  ind^*2*3-2*3 
and  d,.  ■  2  *  3-3*2.  Arc  3  has  zero  flow  and  complementary 
slackness  is  satisfied.  Arc  5  has  f^  »  c,.  and  d^  >  0.  This  fails 
the  test  for  complementary  slackness  and  k,  *  5.  In  this  case, 
flow  will  be  removed  from  arc  5  and  arc  2  and  added  to  arcs  1  and 
4.  As  flow  is  removed,  the  flow  on  arcs  2  and  5  both  reduce  to 

zero  before  the  flow  on  arc  1  or  4  reach  their  upper  bound.  In 
this  case,  arc  5  is  considered  the  limiting  arc  and  k_  *  k.  *  5. 

2j  b 

The  resulting  flows  and  "  updates  are  3hown  in  Figure  5-1 c. 


The  basis  for  Figure  5-1  c  remains  the  same  as  Figure  5-1 b 
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(1,2,4),  with  a  total  coat  of  9.  For  lc_  =*  5 ,  d,  *  2  which  implies 
that  for  aach  unit  change  of  flow  on  arc  5,  the  total  coat  will 
decrease  by  2.  For  a  change  of  2  units,  total  C03t  becomes  13  - 
2(2)  *  9-  Returning  to  step  1  of  the  algorithm  reveals  that  the 
flows  on  the  network  of  Figure  5-1 c  are  optimal. 

Figures  5-2  and  5-3  are  intended  to  3how  the  differences 
between  linear  and  nonlinear  networks  with  respect  to  the  3tepa  of 
the  algorithm.  The  network  of  Figure  5-2  has  the  3ame  form  aa 


Figure  5- 

■1  but 

diffe rent 

arc  parameters.  Arc  1  is  now  a 

nonlinear 

arc  with 

an 

2 

arc  cost  function  of  hf  »  f . .  This  arc  is 

presently  J 

nonbasic 

wi  th 

zero  flow 

as  shown  in  Figure  5-2a. 

The  basis  1 

includes 

arcs 

(2.3,4). 

Step  1  of  the  algorithm  evaluates  the  ' 

nonbasic  area  for  complementary  slackness,  *3+10-3*5  and 
complementary  slackness  is  satisfied.  Arc  1  however,  has  i  =  0  + 
2f^  -  6  =*  -6  since  f1  ”  0.  Thus,  k?  =  1  .  Step  2  of  the  algorithm 
finds  the  arc  which  limits  the  flow  change.  As  flow  is  increased 
on  arc  1  flow  must  be  decreased  on  arcs  2  and  3.  If  the  network  of 
Figure  5-2  were  linear,  it  would  be  profitable  to  add  4  unit3  of 
flow  to  arc  1  and  remove  4  units  from  arcs  2  and  3,  causing  all 
three  of  these  arcs  to  go  to  a  bound,  in  which  case  any  one  of 
these  arcs  could  leave  the  basis.  For  ^  *  -6,  this  would  reduce 
the  total  cost  to  9  3ince  h,  «  0.  The  marginal  cost  for  arc  1  will 
not  remain  zero  as  flow  i3  added  for  the  nonlinear  case. 
Consequently,  it  is  necessary  to  find  the  amount  of  flow  which  when 
added  to  arc  i  will  cause  d^  to  go  to  zero.  In  this  case,  for  f  » 
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3,  hj  -6  and  d1  *  0  +  6  -  6  *  0.  Thus,  the  flow  change  is  limited 

to  3  on  the  nonlinear  arc.  Comparing  this  limiting  flow  with  the 

limiting  flow  for  the  linear  arcs  of  4  and  choosing  the  smallest 

results  in  arc  1  being  the  limiting  arc. 

Figure  5 -2b  shows  the  results  of  this  flow  change  and  the 

updated  node  potentials.  Several  things  should  be  noted  in  Figure 

5-2b.  First,  as  in  the  linear  network,  if  kp  is  the  limiting  arc, 

kT  »  Sc,,.  Thus,  the  basis  of  Figure  5-2b  is  the  same  as  Figure 
L  i 

5-2a.  However,  flow  on  arc  1  (a  nonbasic  nonlinear  arc)  is  not  at 
a  bound.  This  is  allowed  for  nonlinear  networks  but  not  for 
linear  networks. 

Secondly,  when  evaluating  the  for  the  nonbasic  arcs  of 
Figure  5-2b,  d1  *  0.  Thi3  is  due  to  the  previous  iteration  which 
forced  this  result.  Evaluating  arc  5,  yields  d^  »  5  and  the 

network  of  Figure  5-2b  is  optimal.  The  total  cost  is  now  23. 

A  third  feature  to  note  is  that  the  total  cost  did  not 

reduce  to  9  as  it  would  have  if  this  were  a  linear  network.  Thus, 

for  the  nonlinear  network,  the  evaluation  of  d^  is  necessary  to 
test  for  complementary  slackness,  but  it  cannot  be  interpreted  as 
the  change  in  total  cost  per  unit  change  of  flow  in  arc  k„  for 

ii 

nonlinear  networks,  rather  it  is  the  derivative  of  this  change. 

In  the  example  above,  the  node  potentials  did  not  change 
since  the  basis  remained  unchanged  and  contained  only  linear  arcs. 
If  the  nonlinear  arc  had  entered  the  basis  and  since  there  are  no 
other  nonlinear  arcs  in  the  basis,  the  node  potentials  would  be 
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updated  in  the  usual  manner,  that  being  only  to  update  the  node 
potentials  in  the  subtree  rooted  at  the  terminal  node  of  the 
entering  arc.  For  a  more  complicated  network  there  may  be  other 
nonlinear  arcs  in  the  basis.  These  arcs  may  or  may  not  be  in  the 
subtree  rooted  at  the  terminal  node  of  the  entering  arc.  If  they 
are,  the  node  potentials  will  be  updated  routinely.  Basic 
nonlinear  arcs  in  another  part  of  the  network  may  not  experience  a 
flow  change,  however,  if  their  costs  functions  include  an 
interactive  term  relating  to  an  arc  whose  flow  did  change,  then  the 
marginal  cost  associated  with  the  nonlinear  arc  will  change.  This 
requires  that  the  node  potentials  be  updated  for  all  nodes  in  all 
subtrees  rooted  at  the  origin  of  the  nonlinear  arcs  in  the  basis. 

Consider  the  simple  example  of  Figure  5-3.  There  are  two 
nonlinear  arcs  in  this  network,  arcs  1  and  2.  Let  the  cost 
associated  with  arc  1  be  h((fjj)  •  ff  *  .5f^f2  and  fur'ther  suppose 
that  the  flow  on  arc  2  is  zero.  Consequently,  the  node  potentials 
for  nodes  a,b,  and  c  are  (l  ,11,16)  respectively.  Let  k_,  a  2  enter 

th 

the  basis  with  a  flow  of  2  and  suppose  k^  »  arc  3.  The  node 
potentials  for  nodes  d,e  and  f  would  be  updated  in  the  usual  way 
since  they  are  in  the  subtree  rooted  at  the  terminal  node  of  the 
entering  arc.  For  the  linear  problem,  these  are  the  only  node 
potentials  that  would  require  updating.  In  the  nonlinear  case, 
since  the  cost  of  arc  1  is  a  function  of  the  flow  on  both  arcs  i 
and  2,  the  marginal  cost  for  arc  1  will  change. 


Ill 


re  5-3 
Cost  Changes 
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-  1  *  .5(2)  -  2 

Thus,  the  marginal  coat  now  associated  with  arc  1  is  2  instead  of 
t.  Accordingly,  the  node  potentials  for  nodes  a.b  and  c  will 
require  updating  to  values  of  (2,12,17)  respectively.  This  is 
unique  to  the  nonlinear  network  and  must  be  accounted  for  in  the 
algorithms.  The  node  potentials  for  nodes  h  and  i  will  remain 
unchanged  if  arcs  4-  and  5  have  linear  arc  cost  functions. 

The  main  changes  to  the  linear  network  codes  to  allow  use 
of  the  primal  simplex  algorithm  for  the  nonlinear  problem  involve 
three  primary  areas: 

1 .  Determining  the  arc  costs  and  node  potentials 

2.  Calculating  the  effect  of  flow  change 


5.  Finding  the  value  of  flow  change  which  causes  i  to  go 

Zj 

to  zero  for  a  nonlinear  arc  kp. 

These  three  issues  will  be  addressed  in  the  remaining  three 
sections  of  this  chapter. 


*  0m-*. 


■W?  ■■  ■ 
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5.7  Arc  Cost  and  Node  Potential 

The  dual  variable  If  ^  represents  the  marginal  cost  of 
increasing  flow  to  node  i.  For  the  pure  linear  problem: 


(3): 


Tf^i*hk 

Setting  the  node  potential  for  the  slack  node  equal  to  zero  allows 
the  determination  of  the  node  potential  for  all  other  nodes  in  the 
network.  ?or  the  generalized  linear  problem: 


Sq.  (A): 


r.*  77V  hk 


k(  i .  j  )  €.  * 


*k 

For  the  linear  problem  the  marginal  cost  of  flow  through  an 
arc  is  constant.  For  the  nonlinear  network,  the  marginal  cost  is 
not  constant  but  is  a  function  of  flow.  The  marginal  cost  is 
the  partial  derivative  of  the  arc  cost  function  with  respect  to  f^. 
If  nodes  i  and  ,j  happen  to  be  on  a  cycle  a  more  complex  equation  is 
required  for  calculating  the  node  potentials.  Let  M„  » 
0,2 . krt)  represent  the  nodes  of  a  cycle.  The  7f  value  at  node 


1  is: 

Sq.  (“5 ) : 


nr. 


hi 


hk 


k-l 

IT 


j*1 


a  . 
J 


3-1 


where  3  is  the  cycle  gain,  defined  as  the  product  of  the  arc  gains 
on  the  cycle. 


•i..  r  it ;  «****. 


The  numerator  of  equation  5  is  referred  to  as  the  unit  cost  of  the 
cycle.  Once  one  of  the  node  potentials  for  a  node  on  the  cycle  is 
determined,  the  remaining  can  be  calculated  using  equation  4. 

From  equation  2  the  marginal  cost  for  linear  arcs  is  just 
h^  and  for  the  nonlinear  arcs  in  a  quadratic  model  is: 

Sq.  (6): 

hkmV2V» 

Thus,  for  the  nonlinear  problem: 


Sq.  (7): 


nr* 

j  _  fo  r  ( '<(  i .  j )  £•  Mjj  C  Mg 

ak 


5.4  The  Effect  of  Flow  Changes 

When  an  arc  is  to  enter  the  basis,  flow  is  changed  in  the 
basic  arcs.  To  determine  what  these  flow  changes  are,  a  quantity 
is  computed.  ^  represents  the  flow  change  through  node  j,  for  a 
unit  flow  change  in  the  entering  arc. 

The  equations  which  are  presented  here  without 
justification  are  described  in  detail  in  Jensen  and  Barnes  (1980). 
It  is  important  to  present  them  in  this  form  because  in  a  nonlinear 
problem  the  flow  in  basic  arcs  effects  the  node  potentials  and 


hence  the  value  of  i,  .  As  stated  above,  it  is  possible  that  a 

Zi 

flow  change  will  drive  i,  to  zero  before  an  arc  is  driven  out  of 

‘  E 

the  basis  and  before  the  flow  on  arc  k„  reaches  its  bound. 
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Assume  that  an  arc  k^i^.jp)  is  chosen  to  enter  the  basis 
such  that: 

v  Vv%v  °  ■ 

The  quantity  Y.  ia  the  amount  of  flow  change  through  node  j  for 

O 

each  unit  change  of  flow  in  arc  kff.  Depending  upon  the  current 
configuration  of  the  basis,  and  on  the  location  of  node  j  within 
this  configuration,  various  equations  apply  for  calculating  Y.. 

To  facilitate  understanding  of  these  equations,  reference 
is  made  to  Figure  5-4a  and  5-4b.  The  key  to  the  various  equations 
lies  in  the  positioning  of  node  ,j  in  the  basi3  network.  If  the 
basis  contains  a  cycle,  the  basi3  tree  will  be  broken  into  two 
parts  or  3eraitrees  as  shown  in  Figure  5-4.  Assume  that  the  trees 
and  3erai trees  of  Figure  5-4  represent  a  basis  for  3ome  12  node 
network.  In  Figure  5-4a  arcs  A  and  3  cannot  both  be  in  the  basi3 
at  the  same  time  since  there  cannot  be  two  cycles  in  the  same 
component  of  the  basis  network.  However,  for  illustrative 
purposes,  both  are  3hown  here.  The  nodes  indicated  as  i_  and  .j- 

Cl  Cl 

represent  the  origin  and  terminal  nodes  for  the  entering  arc 
respectively.  The  node  notation  Jn  in  the  figure  refers  to  the 
node  in  the  basis  which  corresponds  to  the  equation  numbers  below: 


_ ti 


(b) 

Figure  5-* 
Basis  Trees 
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Sq.  (3):  Wien  node  j  is  on  the  basis  path  to  node  i^  but  not  on  a 
path  to  jg  nor  on  a  cycle: 


*>• 


TTS 

k£P..  K 
J1S 

where  P..  i3  the  directed  path  from  node  j  to  node  i_  defined  by 
the  basis.  Por  a  basis  network  this  path  is  always  unique.  Node 
J8  in  Figure  5-4-a  is  an  example  of  a  node  for  which  equation  (3)  is 
appropriate . 

Sq.  (9):  If  node  j  is  on  the  basis  path  to  node  ig,  not  on  a  path 
bo  jp,  but  lies  on  a  cycle: 


(B-1 ) 


where  3  is  the  cycle  gain.  To  illustrate  equation  9  in  Figure 
5-4a,  remove  arc  8  from  the  basis  and  add  arc  A.  Now,  J9  is  a  node 
corresponding  to  equation  (9).  For  equations  ( 1  0 )  and  ( 1  1  )  refer 
to  the  original  basi3. 

Sq.  (l 0) :  For  nodes  such  that  there  exists  a  directed  path  from 

node  j  to  node  (using  basic  arcs)  but  not  to  i_  and  node  j  is 

D  u 


I 

j 

I 

| 


not  on  a  cycle: 


5 

J 


TT’-< 


where  P..  is  the  set  of  area  on  the  oath  from  node  ,j  to  node  j_. 

=< 

Eq.  (11):  Por  nodes  such  that  there  is  a  directed  path  from  node  j 
to  node  .jj,  (using  basic  arcs)  but  not  to  i^  and  node  j  is  on  a 


cycle: 


(B~t  ) 


k£p: . 


Sq.  (12):  If  there  exists  a  directed  oath  from  node  .j  to  both  i„ 

and  end  j  is  not  on  a  cycle: 


TT 

lcfcP,. 


This  is  equal  to  equation  (9)  plus  equation  (10). 
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3q.  (13):  If  there  exists  a  directed  path  from  node  .j  to  both  iff 

and  j„,  and  j  happens  to  be  on  a  cycle: 

Si 


TT 

k£P.. 

•Hr 


IT  V 


!c  €.P  .  . 

J  J51 


This  is  equal  to  equation  (9)  plus  equation  (10). 


Nodes  that  do  not  lie  on  a  directed  oath  to  either  i„  or 

-  £j 


,j?  have: 


This  i3  represented  by  nodes  X  and  Y  in  Figures  5-4a  and  5-4-b. 

Note  that  a  positive  value  of  Y.  indicates  the  flow  in  basic  arc 
Vc( i . j )  will  increase,  and  a  negative  value  indicates  that  thi3  flow 
will  decrease. 

If  arc  'ic(i.j)  is  a  member  of  the  basis  the  flow  change  in 

arc  k  per  unit  increase  in  f,  is: 

3 

3q.  (15): 


For  convenience  define  the  quantity: 


a, 

k 


k(i,j)£  M 

B 
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The  maximum  flow  increase  for  arc  kv  is  that  value  that  will  cause 
a  basic  arc  (or  arc  ~<v)  to  go  to  a  bound.  This  maximum  is: 


Min[(ck  -  fk)  gk  k(i.j)£M3,  gk  >  0]  . 

Min  [-fk  gk  k(i..j)€M3,  gk  <  o] 

Thi3  equation  is  used  if  d,  <0  and  f  <  c.  .  The  maximum  flow 

*3  *3 


decrease  in 
3q.  (17): 


Min  Qk  ^k  ^3-  ?)c  >  > 


“ 

Q^ck  "  Mp  «k 

k(i,.j)  €  Mg,  gk  < 

This  equation  is  used  if  d,  >0  and  f,  >0.  Equations  (IS)  and 

*3  *3 

(17)  are  U3ed  to  determine  the  flow  change  which  will  drive  an  arc 
flow  to  a  bound.  For  the  nonlinear  problem  it  is  also  necessary  to 
determine  the  flow  change  that  will  cause  d,  to  go  to  zero.  Thi3 

*3 


I 

\ 


I 


is  the  subject  of  the  next  section. 
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5.5  Flow  Change  Which  Drives  d,  to  Zero 


For  the  nonlinear  problem  the  value  of  4.  is  determined 


by: 

Eq.  (IS) 


^i  *  hk  '  ak  "A \ 
LS  E  E  JE 


Assume  f,  <  c.  and  d,  <  0  indicating  a  nonoptimal  condition. 
<E  o  ^E 

The  results  of  the  last  section  indicate  the  flow  change  that 
causes  a  basic  arc  (or  arc  !<„)  to  go  to  a  bound.  Since  7T.  ,  h' 

i  i-ri 

Zj  Zj 

and  IT  .  may  all  change  as  flow  changes  in  a  nonlinear  problem,  it 
JE 

is  possible  that  a  smaller  flow  change  than  that  determined  in  5-4 

will  cause  d.  to  go  to  zero.  If  this  were  to  hapoen,  arc  <_  would 
SE  1 

no  longer  violate  complementary  slaclcness.  Thi3  section  derives 

the  value  of  flow  change  in  arc  which  drives  d,_  to  zero. 

An  alternative  form  of  d,  i3  useful  here.  The  value  d, 

i3  the  marginal  cost  change  with  respect  to  a  flow  change  in  ar 
k  Since  the  total  cost  changes  only  with  the  flow  changes  on 

basic  arcs  occasioned  by  the  flow  change  on  arc  the  value  of 


d,  can  be  written: 


Sq.  (19): 


k(i.j)  £  \ 


Hote  that  g^  is  the  marginal  flow  change  in  arc  k  with  respect  to  a 


flow  change  in  arc  k,  and  h(  is  the  marginal  cost.  Since 


o  lU 


122 


h,^=h,<+2QkF?(,  separating  the  linear  and  nonlinear  terms  yields: 
3q.  (20): 


i, 

<, 


h!<gk  ^^k 


3  3 


(2QkF^)gk 


This  sum  is  over  all  the  arcs  in  the  basis  since  the  value 


g  is  zero  except  for  areas  on  paths  which  are  effected  by  a  flow 

change  in  arc  In  the  linear  network,  d,  is  not  effected  as 

flow  changes  in  arc  3ince  h^  is  independent  of. flow.  However, 

if  a  nonlinear  arc  k(i,j)  has  g^O,  changing  flow  effects  h,(,  and 

hence  d.  .  What  one  would  like  to  be  able  to  do  is  to  compute  the 
*3 

value  of  flow  change  in  the  entering  arc  at  which  d.  goes  to  zero. 
Define: 

3q  •  (21): 


Ad 


where  d,  was  the  value  before  the  flow  change. 

.Igain  separating  the  linear  3nd  nonlinear  terms: 
3q.  (22): 


d,' 

k_ 


3  kcM, 


hkgk  *  ? 


\r  *\C  *  - 

*3  *3 
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In  thi3  equation,  is  the  vector  of  flows  in  the  nonlinear  arcs 
modified  as  flow  in  arc  k_  is  changed. 

Let: 

3q  •  (23): 

Ps-V  Afjj 


I 


where  F^  is  the  original  flow  and  ^  F^  is  the  incremental  flow  in 
the  nonlinear  arcs. 

Sq.  (24): 


K  ^  h,  g.  *  h'  g.  *  (29,  F„)g. 

*5  k£«B  *  *  *8  *3  k£MN  *  *  * 


(2Q,^F„)g 
kfcM^  <  *  < 


The  first  three  terms  in  this  expression  comprise  d,  <md  the  las* 
term  is  ^d . 


Solving  for  the  flow  change  that  makes  1).  *0: 


Sq.  (25): 


&d=-d. 


Vs£r  (2tW*« 


k6M„ 


Define  for  the  nonlinear  arcs  the  vector  3„  where: 

N 


Sq.  (26): 


for  k  £ 


Since  9  is  symmetric,  equation  25  can  be  rewritten: 


Sq.  (27): 
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From  the  previous  discussion: 

Bq.  (23): 

^P!f=5Nfk_ 

E 


Substituting  for  F^  j.n  (27)  yields: 
Eq.  (29): 


Solving  for  the  flow  ohange  that  causes  d'  to  be  zero: 
Sq.  (SO): 


9  h 

where  A?.  is  the  change  in  flow  f,  . 

kE  kE 

Note  that  3ince  d,  is  negative  and  Q  is  positive  definite, 
*E 

the  value  of  /\f,  will  be  oositive.  For  the  case  where  d.  is 
E  '  KS 

positive,  a  negative  value  for  f,  results.  This  is  appropriate 

S 

since  it  indicates  that  the  flow  on  arc  k-  must  be  reduced  to  drive 


d,,  to  zero. 

*E 

The  result  of  the  above  formulation  is  a  value  for  flow 

change  in  the  nonlinear  arc  k7  which  causes  d,  to  be  zero.  This 

J  *E 

limiting  value  is  compared  with  the  flow  change  determined  by 
equation  (16)  or  (17).  The  minimum  of  these  two  values  will 

determine  the  amount  of  flow  change  to  apply  to  arc  kP.  If  the 
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limiting  arc  is  determined  by  equation  (16)  or  (17)  arc  kp  will 
enter  the  basis  and  the  limiting  arc  will  leave.  If  the  limit  i3 
obtained  by  equation  ("30),  the  flow  will  be  changed,  but  the  basis 
will  remain  unchanged.  In  this  case,  the  nonbasic  nonlinear  arc 
will  have  flow  between  its  bounds. 

The  above  represents  a  general  development  for  solving  the 
nonlinear  nonseparable  quadratic  network. 

This  presentation  justifies  the  inclusion  of  these 
nonlinear  quadratic  arcs  in  the  network  in  their  nonseparable  form. 
To  implement  this  new  theory  into  the  linear  codes  of  Jensen  and 
Barnes  (i960),  six  new  subroutines  were  created  and  several 
existing  subroutines  required  some  modification.  These  new  codes 
were  assembled  in  a  package  called  MONLING  and  tested  in  the 
presence  of  positive  and  negative  gains.  These  subroutines  are 
flow  charted  in  Part  III  of  the  Appendix. 
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CHAPTER  VI 

6.  Estimating  the  Benefit  Function 

5.1  General 

The  preceding  sections  have  shown  that  by  exploiting  the 
special  structure  of  networks  in  a  dynamic  programming  approach  an 
expression  can  be  derived  which  represents  the  future  expected 
value  of  stored  water.  This  was  done  by  utilizing  a  single  period 
network  representation  of  a  water  resources  system.  Through  the 
selection  of  discretized  initial  levels  and  the  incorporation  of 
stochastic  inflows  to  these  initial  levels,  the  total  water 
available  to  the  period  was  determined.  Subsequently,  through 
successive  random  draws  from  the  distribution  of  the  stochastic 
inflows,  enough  data  was  assumed  to  be  generated  to  allow  a 
reasonably  good  approximation  of  a  benefit  function  using  a  least 
squares  fit  of  the  observed  network  optimal  solutions. 

The  water  resources  problem  solution  procedure  requires  two 
types  of  information  in  order  to  produce  meaningful  and  useful 
results.  The  first  type  relates  to  the  specific  network 

parameters.  These  include  all  of  the  arc  parameters,  the  river  and 
reservoir  data,  demands,  runoff  distribution,  etc.  Naturally,  this 
information  is  problem  specific  and  would  most  generally  be 
supplied  by  the  user.  A  second  class  of  information,  however,  is 
both  model  specific  and  user  oriented  and  relates  to  the  accuracy 
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and  credibility  of  the  model.  Thi3  class  requires  the  answer  to 
several  questions. 

1 .  How  many  discretizations  of  reservoir  contents  should  be 
used  and  from  the  total  number  of  possible  level 
combinations  generated,  which  should  be  used  for  an  optimal' 
design  for  this  experiment? 

2.  What  method  of  least  squares  regression  should  be  used 
to  fit  the  data  obtained? 

3.  How  many  replications  from  the  assumed  ’mown  inflow 
distribution  should  be  used  for  each  level  combination? 

4.  How  can  convexity  be  assured? 

After  addressing  these  four  questions,  and  having  observed  some  of 
the  results,  additional  concerns  require  attention.  These  include: 

5.  Examining  the  variance  covariance  matrix  to  measure  to 
validity  of  the  least  squares  regression 

6.  Consideration  of  weighted  least  squares 

7.  Problems  associated  with  weighted  least  squares 
Finally,  based  upon  all  of  the  above  analysis: 

3.  Determination  of  the  data  requirements  for  estimating 
the  benefit  function. 

These  questions  and  related  issues  will  be  addressed  in  the 
next  3  sections. 

6.2  Design  of  the  Experiment 


The  process  of  estimating  the  benefit  function  requires  the 
discretization  of  the  water  levels  in  the  system  reservoirs. 
Consider  a  three  reservoir  system  where  each  reservoir  has  a 
capacity  of  25  units  of  water.  Each  unit  will  represent  3orae 
specified  number  of  acre  feet.  The  hypothetical  reservoir  system 
of  Figure  4-1  will  be  used  throughout  this  discussion.  This 
hypothetical  example  is  further  exercised  in  detail  in  Chapter 
Uote,  there  is  no  requirement  for  these  reservoirs  to  be  equal  in 
size  and  the  reservoir  sizes  are  input  parameters  to  the  model,  set 
by  the  user  as  dictated  by  his  specific  system. 

Reservoirs  are  typically  divided  into  zones  according  to 
their  functional  purpose.  Linsley  and  Franzini  0964)  define  these 
zones  as  shown  in  Figure  6-1 .  For  a  very  low  level  of  water  there 
is  a  minimum  volume  referred  to  as  the  dead  storage  pool.  This  is 
the  level  below  the  sluiceway  where  water  cannot  be  released  from 
the  reservoir  for  other  uses.  At  the  very  top  of  the  reservoir,  a 
flood  pool  is  normally  reserved.  The  purpose  of  this  pool  is  to 
allow  the  temporary  storage  of  water  during  heavy  inflow 
conditions,  thus  preventing  high  water  flood  conditions  downstream. 
During  normal  operations,  This  flood  pool  will  usually  be  empty. 
Between  these  two  pools  is  the  useful  storage  volume.  It  is  this 
portion  of  thereservoir  that  is  of  primary  interest  to  this  study. 
Thi3  useful  storage  is  the  portion  of  the  reservoir  that  will  be 
discretized  for  the  dynamic  programming  approach.  The  percent  of 
reservoir  volume  in  the  dead  storage  pool  and  in  the  flood  pool  is 
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dependent  upon  the  specific  reservoir.  Within  this  useful  storage 
range,  members  of  the  Texas  Water  Development  Board  indicate  that  a 
reservoir  holding  only  50^  capacity  could  be  considered  as  having  a 
severe  shortage  (depending  upon  the  size  and  use  of  the  reservoir). 
Thus,  although  the  total  useful  storage  pool  is  available,  the  user 
may  be  concerned  about  only  a  portion  of  it.  Thus,  the  range  to  be 
discretized  will  be  left  to  the  user  and  he  has  complete 
flexibility  (within  the  usable  storage  pool),  in  determining  this 
range  3ince  the  maximum  and  minimum  water  levels  are  input 
parameters  to  the  model  as  a  percent  of  reservoir  capacity. 
Pegardless  of  the  range  it  is  necessary  to  determine  the  number  of 
levels  or  discretizations  to  use  to  produce  a  meaningful  set  of 
data,  and  finally,  what  subset  of  the  total  set  of  level 
combinations  can  be  used  to  achieve  an  optimal  design? 

The  significance  of  the  number  of  discretizations  lies  in 
its  effect  on  the  accuracy  of  the  quadratic  fit  and  on  the 
computational  time  required.  Ideally,  one  should  use  as  few 
discretizations  as  possible  which  yield  the  best  "or  near  best”  set 
of  data  for  a  quadratic  fit.  The  total  number  of  possible  level 
combinations  is  equal  to  the  number  of  levels  (NN)  raised  to  the 
number  of  reservoirs  (R)  power.  Jfotationally ,  the  number  of  level 
combinations  is: 

rmR 


regardless  of  the  range. 

For  a  three  reservoir  problem  with  three  levels  selected 
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*  27  combinations  of  initial  reservoir 

five  levels  there  are  5^  *  125  total 
this  to  a  larger  problem,  say  six 
reservoirs  and  five  levels  there  would  be  1 5625  combinations. 

As  noted  in  the  dynamic  programming  algorithm  of  Chapter  4, 
random  inflows  are  added  to  the  reservoir  levels  at  the  beginning 
of  the  period  to  yield  total  water  available.  Network  solutions  to 
these  problems  yield  the  desired  observations  which  are  used  to  fit 
a  quadratic.  Naturally,  one  does  not  need  15625  data  points  to  fit 
a  quadratic  having  29  coefficients  (which  is'  the  case  for  a  six 
reservoir  problem) . 

Assume  for  now  the  three  reservoir  case  and  three  levels 
of  discretization  for  each  reservoir.  The  problem  then  is  to  select 
from  the  total  set  of  level  combinations  an  optimal  3et  of  design 
points.  Using  this  terminology,  a  design  point  relates  to  a  level 
combination  Figure  6-2  shows  a  three  level  cubic  representation 
for  the  three  reservoir  problem.  Note,  there  are  27  distinct 
points  on  and  within  this  cube  represented  by  the  intersection  of 
lines  (ignore  the  highlighted  points  for  now).  Each  of  these 
points  represents  a  level  combination  or  a  design  point.  The  goal 
is  to  select  some  subset  from  these  2^  design  points  which  will 
yield  an  optimal  design.  3ox  and  Draper  (1971)  and  Mitchel  ( 1974a, 
1974b)  have  defined  an  optimal  set  of  design  points  for  first  and 
second  order  models.  Their  work  is  specifically  for  least  squares 
approximations  where  the  underlying  error  has  a  normal 


for  each,  there  are  3^ 
levels  considered.  For 
combinations.  Extending 


< 


* 


I 


Hu  *  ■*  vJ.  -  kOn,, 


distribution. 


Two  of  the  basic  approaches  to  choosing  an  n-point 
experimental  design  are  (i)  to  set  down  a  simple  factorial  or 
fractional  factorial  design  in  the  factors  being  studied,  or  (ii) 
to  choose  a  design  based  on  the  well  known  det(X'X)  criterion. 
Box  and  Draper  (1971)  indicate  that  the  first  method  is  much  more 
simplistic  and  that  the  second  is  their  preferred  approach. 

Details  of  their  preferred  approach  are  presented  below  and 
have  been  utilized  in  the  computer  program. 

3efore  presenting  their  approach,  a  few  new  terms  need  to 
be  discussed.  A  design  point  has  already  been  defined  as  being 
equivalent  to  a  level  combination.  A  design  point  has  dimension 
’ ixH)  for  an  R  reservoir  problem.  The  design  matrix  was  defined  in 
Chapter  1  as  being  the  matrix  of  independent  variables  to  be  used 
for  the  regression.  A  full  quadratic  has  been  selected  as  the 
model  to  be  fit  by  the  regression.  For  the  three  reservoir,  three 
discretization  problem,  this  yields  a  design  matrix  of  dimension 
,?"’xiO)  if  all  design  points  are  used.  To  use  the  notation  of  Box 
and  Draper  and  to  be  consistent  with  standard  notation  for 
regression  analysis,  this  design  matrix  will  be  referred  to  as  the 
X  matrix. 

In  their  preferred  approach,  the  criterion  for  designing 
experiments  is  based  on  maximizing  the  determinant  of  (X'X) 
indicated  here  as  det(T'X),  (X'  is  the  transpose  of  the  X  matrix). 
TJse  of  this  criterion  dates  back  to  Smith  (1913).  The  criterion 
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has  many  appealing  properties  and  its  use  has  been  .justified  in  a 
number  of  ways.  Box  and  Lucas  (1959)  indicated  that  its  use  leads 
to  a  confidence  region  for  the  parameter  estimates  of  smallest 
hypervolume  in  parameter  space.  Kiefer  ( 1 961 )  showed  that  a  design 
which  maximizes  det(X'X)  also  minimizes  the  maximum  variance  of  any 
predicted  value  (obtained  by  using  the  regression  function)  over 
the  experimental  3pace.  Further  properties  of  thi3  criterion  are 
that  it  minimizes  the  generalized  variance  of  the  parameter 
estimates,  and  that  the  design  obtained  is  invariant  to  changes  of 
scale  of  the  parameters.  This  is  an  important  property  not  shared, 
for  example,  by  the  criterion:  minimize  trace  (T'X)”',  (i.e.,  min 
the  average  variance  of. the  parameter  estimates). 

Box  and  Draper  (1971)  considered  several  discretization 
techniques  and  found  that  the  best  design  for  a  quadratic  fit  was 
to  use  three  levels  consisting  of  both  end  points  and  the  midpoint 
of  the  selected  range.  These  design  points  correspond  to  the  27 
points  of  Figure  6-2.  Thus,  there  are  27  design  points  from  which 
to  choose. 

The  choice  of  three  levels  of  discretization  runs  contrary 
to  the  typical  dynamic  programming  approach.  Klemes  and  Doran 
(1977)  specify  5-10  levels  in  their  divided  interval  technique, 
while  others  3ay  as  many  as  30  levels  may  be  needed.  However, 
these  levels  are  not  being  used  in  the  typical  dynamic  programming 
manner.  The  requirement  is  to  select  enough  levels  to  allow  a  good 
approximation  of  a  quadratic  benefit  function  and  three  points  are 
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sufficient  to  fit  a  quadratic-  In  practice,  once  an  expression 
which  reflects  the  future  value  of  water  has  been  derived,  the  user 
can  apply  the  actual  water  levels  to  this  function  resulting  in  an 
operating  policy  representative  of  the  current  situation. 

What  a  ust  be  done  next  is  to  select  n  of  these  design 
points  and  apply  the  det(X'X)  criterion,  (n  in  this  case  can  vary 

from  10  to  27).  The  design  matrix  for  the  full  quadratic  would  now 

27 

have  dimension  (nxlO).  For  any  given  n,  there  are  C£  subsets  to 

consider,  and  for  each  n,  there  will  exists  a  "best"  design  set  in 

terms  of  maximizing  det(X'X).  For  example,  if  n  =  20  there  are 
27 

CpQ  *  388,030  different  subsets  of  20  design  points  to  evaluate. 

Oalculating  the  det(X'X)  for  all  388,030  of  these  subsets  will 
result  in  one  of  them  having  a  determinant  that  is  greater  than  or 
equal  to  all  others.  The  subset  with  the  maximum  det(X'X) 

represents  the  best  set  of  design  points  to  use,  given  that  20  are 
to  be  used.  Now,  plotting  the  optimal  det(X'X)  obtained  for  each  n 
yields  a  form  similar  to  the  one  shown  in  Figure  5-3.  Through 
this  process.  Box  and  Draper  were  able  to  derive  an  optimal  design 
for  a  quadratic  which  has  a  general  structure  and  can  be  applied  to 
problems  of  three  or  more  factors. 

The  resulting  optimal  design  is  referred  to  as  the  "cube 
plus  3tar"  design.  Mitchel  (1974a,  19'1-lb)  refers  to  it  as  ”D 
optimal"  experimental  design  (the  D  referring  to  the  determinant  of 
the  X’X  matrix).  This  "cube  plus  star”  nomenclature  essentially 
describes  the  design.  In  Figure  5-2,  the  "cube"  refers  to  the 
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eight  vertices  of  the  cube.  The  "star"  i3  defined  as  the  center 
point  on  each  of  the  3ix  surfaces  of  the  cube  and  one  point  exactly 
in  the  center  of  the  cube.  This  results  in  a  total  of  15  design 
points  of  a  possible  27  points  in  the  three  factor  case.  These  15 
points  are  the  highlighted  points  of  Figure  6-2. 

In  general,  if  there  are  R  reservoirs  or  dimensions  to  the 
problem,  there  will  be: 

2R  *  2R  +  1 

design  points.  This  means  that  for  the  six  reservoir  case 
mentioned  before,  rather  than  using  five  discretizations  for  a 
total  of  15625  design  points,  only  three  discrete  levels  and  77 
total  design  points  are  required. 

One  possible  disadvantage  of  the  det(X'X)  criterion  is  that 
it  is  a  “variance  criterion"  and  effectively  assumes  that  the  model 
considered  is  the  true  model.  Box  and  Draper  ( 1 971 )  point  out 
however  that  in  situations  where  the  design  is  physically 
restricted  to  a  cuboidal  region  of  interest,  the  difference  between 
the  spread  of  the  design  points  for  the  best  all-bias  design  and 
for  the  best  all-variance  design  is  minimal.  Thu3,  they  conclude, 
that  the  det(X’X)  criterion  appears  to  be  not  unrealistic  either 
when  the  model  is  correct  or  when  the  design  is  restricted  to  the 
region  of  interest,  or  both. 

The  det(X'X)  criterion  applies  to  normal  least  squares 
problems  where  the  variance  about  each  design  point  is  constant  and 


there 


13  no  oovsrisnce  present. 


Thus,  the  variance  covariance 
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matrix  has  the  standard  form: 

V 

0 
0 

The  variance  covariance  matrix  for  the  three  reservoir 
problem  will  have  dimension  (15x15)  where  each  row  represents  one 
of  the  15  design  points.  The  columns  have  the  3ame  meaning.  The 
values  on  the  main  diagonal  (<T~)  represent  the  variance  about  the 
mean  of  the  calculated  response  for  each  design  point.  This 
variance  is  determined  by  taking  a  number  of  replications  at  each 
design  point  and  calculating  the  variance  using  standard 
statistical  methods.  This  variance  is  considered  to  be  the  same 
for  all  design  points.  This  is  the  definition  for  homoscedasticity 
of  data,  Clark  and  Schk3de  (1974).  All  variance  covariance 
matrices  are  symmetrical.  This  particular  one  has  all  off  diagonal 
elements  equal  to  zero.  This  indicate^  that  there  is  no  covariance 
present.  This  is  necessary  for  the  use  of  normal  least  squares 
regression  which  will  be  presented  next. 
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6.3  Normal  Least  Squares 

In  the  analysis  of  variance,  for  the  normal  least  squares 

method : 

b  -  (rxr'x'Y 

where: 

X  is  the  design  matrix  (independent  variables)  X ' *X 
transpose 

Y  i3  the  vector  of  dependent  variables  (mean  of  the 
replications) 

b  is  the  vector  of  estimates  of  the  coefficients 
The  variance  of  the  coefficients  is  found  by: 

Sq  (t): 

7(b)  »  (X'X)"1  (X'VXHX’X)'1 

where  V  is  the  variance  covariance  matrix.  This  reduces  to: 

V(b)  =«  (X'X)"1  O'2 

p 

when  (T  *  is  icnown  and  7  is  considered  to  be  the  Identity  matrix. 
For  independent  observations,  this  assumes  that  the  variance 
covariance  matrix  of  the  design  matrix  is  of  the  form: 


Here  as  before,  the  diagonal  element  (i,i)  represents  the  variance 
of  the  replications  about  the  mean  of  the  ith  design  point.  For 
the  water  resources  problem,  is  not  known.  Consequently,  an 
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estimate  of  <T  must  be  derived  by  performing  enough  replications 
at  each  design  point  to  get  a  reliable  estimate  of  which  in 
thi3  case  is  synonymous  with  the  V  matrix.  Choosing  the  number  of 
replications  is  the  subject  of  the  next  section.  (Jiven  that  ’/  has 
been  derived  of  the  form  shown  above,  the  variances  of  the 
estimated  coefficients  can  be  calculated  using  equation  1 . 


o.4  Replications 

A  means  for  determining  the  number  of  draws  for  deriving  a 
o 

good  estimate  of  “  i3  to  select  some  acceptable  relative  error 
(alpha)  which  represents  the  ratio  of  the  standard  error  of  the 


estimate  of  the  variance  to  the  true  variance. 


It  is  known  that: 


^if5X(W)^Var(^’(r))  *  2r 

(T 2 

p 

where  1c  is  the  number  of  draws,  X  ~  is  the  notation  for  the  chi 
square  distribution  and  Y  is  the  iegrees  of  freedom  for  the  chi 


square  distribution. 


(ic-1  )2  Var(s2)  -  2(k-1  ) 
Var( 3?)  -(2  <J"f)  /  (ic-1  ) 


3E(s2)  -JX  <j  2  £  alpha 


1 


3E  is  the  standard  error  of  the  estimate  of  the  variance,  s^. 
Thus  the  relative  error  of  the  estimate  to  the  true  error  13: 


3E(s2)  n 


For  fixed  values  of  alpha  k  can  be  determined.  A  few 
values  are  shown  in  Table  6-1 . 


Table  6-1 


Number  of  Draws  Required  to  Meet  Specified  Accuracies 
alpha  k 


If  the  process  described  by  the  reservoir  system  were  relatively 
stable  or  unchanging  from  period  to  period,  V  could  be  estimated  to 
the  desired  accuracy  one  time  using  a  large  number  of  draws.  Then 
through  the  use  of  ordinary  least  squares  equations,  the  variance 
of  the  coefficients  of  the  benefit  function  for  each  period  could 
be  determined.  For  this  problem,  many  things  change  from  period  to 
period,  specifically,  the  inflows  and  arc  parameters  (including  the 
3  matrix)  for  the  nonlinear  arcs.  These  changes  require  that  the  V 
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matrix  be  reevaluated  for  each  period  since  each  period  is  in 
effect  a  new  problem.  This  would  suggest  using  as  few  draws  as 
possible  to  keep  the  model  computationally  feasible.  30  draws  will 
be  used  for  the  example  problems  and  this  is  considered  adequate 
for  the  water  resources  problem.  If  the  user  desires  greater 
accuracy,  naturally  he  can  increase  the  number  of  draws  at  the 


expense  of  time. 
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6.5  Maintaining  Concavity  of  Benefit  Functions 

A  requirement  of  the  optimization  procedure  is  that  each 
benefit  function  derived  by  the  quadratic  fit  and  subsequently  used 
in  the  iynamic  programming  procedure  must  be  a  concave  function. 
Thus  for  each  t,  BF( t)  must  be  concave.  This  requirement  comes 
from  a  limitation  of  the  network  flow  optimization  algorithm  of 
Chapter  5.  Specifically  the  algorithm  only  works  to'  find  a  flow 

solution  which  minimizes  total  cost  if  the  total  cost  function  is 
convex.  Since  the  negative  of  the  benefit  function  is  part  of  the 
cost  function,  this  requires  that  the  benefit  function  be  concave. 

This  limitation  of  nonlinear  minimization  algorithms  to 
convex  objective  functions  is  not  uncommon.  The  presence  of 
concave  portions  of  the  objective  function  may  result  in  local 
minimums.  Algorithms  to  handle  the  more  general  problem  are 
usually  much  more  complex  and  require  more  computation  time  for 
obtaining  a  solution.  At  any  rate  the  network  flow  algorithm  can 
handle  only  convex  cost  functions  (or  concave  benefit  functions). 

The  requirement  of  concave  benefit  functons  does  not  3eem 
to  impose  serious  practical  limitations  for  the  water  resources 
problem.  It  is  clear  that  the  marginal  value  of  water  stored  for 
the  future  should  be  declining  with  the  amount  of  water  stored. 
This  can  be  proved  to  be  true  for  the  models  used  for  this 
research.  However  it  is  quite  possible  to  happen  upon  a 


non-concave  quadratic  fit  if  proper  precautions  are  not  taken. 
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Sven  if  the  underlying  model  has  a  concave  benefit,  if  independent 
random  observations  are  taken  at  each  design  point  statistical 
error  may  result  in  a  convex  fit.  This  is  especially  likely  if  the 
model  is  nearly  linear.  One  thing  that  could  be  done  is  to  check 
the  quadratic  matrix  for  negative  definiteness  after  each  benefit 
function  is  derived.  However,  in  the  event  that  the  benefit 
function  is  not  concave,  an  alteration  of  the  fit  would  be  required 
in  order  to  continue.  It  is  not  clear  how  to  perform  this 
alteration  in  the  general  case.  Instead  of  resorting  to  this 
manipulative  alteration,  a  random  sampling  procedure  will  be  used 
which  will  insure  a  concave  form. 

is  indicated  before,  a  series  of  random  numbers  which  are 
used  to  derive  the  stochastic  inflows  is  selected.  If  a  different 
random  number  is  selected  for  every  design  point,  there  is  no 
assurance  that  the  final  fit  will  be  a  concave  quadratic, 
regardless  of  the  number  of  replications. 

To  assure  concavity  the  same  set  of  random  numbers  will  be 
applied  to  every  design  point.  Thus,  if  there  are  to  be  k 
replications  for  each  design  point,  there  will  be  a  total  of  k 
random  numbers.  These  same  k  random  numbers  will  be  used  to 
generate  the  k  replications  for  all  levels. 

To  illustrate  this  idea,  refer  back  to  Figure  4-4..  In  this 
figure,  there  are  design  points.  The  response  matrix  is  (L^xk) 
where  k  i3  the  number  of  draws.  Note  that  the  w^  draw  is  applied 
to  all  1  levels  of  the  design.  Next,  the  Y\  ^  are  averaged  over  the 
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k  replications  and  these  mean  responses  are  then  used  to  fit  the 
quadratic.  Thus,  it  will  be  shown  that  by  using  the  mean  response 
values,  obtained  through  the  application  of  a  constant  3et  of 
random  numbers,  concavity  will  be  assured. 

Consider  a  particular  design  point  i  defined  by  the 

reservoir  contents  S./,.  ,  A  random  draw  w  determines  an  inflow 

U  t-1  ) 

vector  I  .  Since  the  inflows  are  assumed  to  occur  at  the 
w 

reservoirs,  the  total  water  inputs  to  the  system  for  design  point  i 
and  random  draw  w  is: 

3i(t-l)  rw 

These  inflows  appear  in  the  network  model  as  positive  external 

* 

flows  for  the  reservoir  nodes.  The  optimum  flows  are  determined 
for  the  network  model  and  the  minimum  cost  is  the  response  Y. 

i.  *  w 

The  network  problem  has  a  convex  objective  function  since  the 

problem  is  linear  except  for  the  convex  function  -BF(t). 

Let  Y.(l  )  be  the  value  of  the  minimum  cost  solution  of  the 
i  w 

network  model  as  a  function  of  the  vector  I  .  It  is  well  known 

w 

that  the  minimum  value  of  a  convex  objective  expressed  as  a 

function  of  the  right  hand  sides  of  the  constraints  is  also  convex. 

Thus,  Y.(l  )  is  a  convex  function  with  respect  to  I  .  The 
i  w  w 

procedure  used  to  obtain  BF(t-l),  samples  from  the  distribution  of 

I  to  obtain  k  distinct  values.  The  k  response  values  of  Y.  thus 
w  r  1 ,  w 

obtained  must  fall  on  the  function  Y^(l^).  By  following  the  above 
procedure  for  the  water  resources  problem,  the  possibility  of 
generating  response  data  that  was  nonconvex  was  removed. 
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6.6  Covariance  Matrix  Derivation  and  Analysis 

Because  of  the  requirement  to  use  a  constant  random  number 
set  for  all  design  points,  a  high  degree  of  correlation  between 
the  design  points  has  been  induced  into  the  model.  The  fact  that 
this  correlation  is  present  requires  the  consideration  of  using  a 
weighted  least  squares  regression  analysis. 

Whether  weighted  least  squares  or  ordinary  least  squares  is 
used,  it  is  required  that  a  very  good  estimate  of  the  variance 
covariance  matrix  (V)  be  derived.  This  is  necessary  3ince  the 
variance  associated  with  the  problem  is  not  known  apriori.  Because 
of  the  possibility  of  correlation  being  present  it  is  necessary  to 
derive  this  V  matrix  for  use  with  both  the  ordinary  and  weighted 
least  squares  analysis. 

The  set  up  for  weighted  least  squares  analysis  when  all 
design  points  have  common  random  inputs  goes  as  follows: 

Let : 

Y  *  response  at  the  i  design  point  when  the  random 
input  is  the  wth  random  sample  from  a  given  distribution 

Xi  *  'xi1'  xi2’"*’xiu^  "  u“V9C'tor  settings  of  the 

independent  variables  at  the  itfl  design  point,  1  <_  i  <_  n. 
where  n  is  the  number  of  design  points  used  in  the  n-point  optimal 
design. 

The  assumed  model  where  3  is  the  assumed  coefficient  is: 


Y. 

r.w 


where: 


8.  vs(0,  IT2),  1£  i  £  n 

x ,  w  i 

2 

and  where  depends  on  l. 

.Vow  by  supplying  the  same  random  input  w  to  all  n  design  points, 


the  overall  model  in  this  situation  becomes: 


V, 


1  X  X 

,,2  >  S’ Oov(£.  ,  £.  ) 

K  Ofrr  “frr  i.w  3 .« 


?Iow  w  ¥  m  ^  w  and  m  are  independent 

Cov(£.  £..  )  =  E(  £ .  £.  ) 

l.w’  J  ,m  l.w  j , m 

-  S(  t.  )  S(  £,  m) 

l.w  J  .ffl 

The  last  two  terns  here  are  equal  to  zero  by  the  assumption  of  the 
nodel . 


This  leaves: 


3(  £  )  3(  £ i  _)  =  0 

i.w  J .  m 


since  v  and  m  are  indenendent . 


On  the  other  hand,  for  w  =  m: 


Cov(  £  £  )  *  Cov(  i  l  ) 

J  *  ni  l .3  *  ^ 


Cov(  I  .  ,  I  . )  -  I  V°  . 

1  J  X2^  1J 


1  , ,  t.O  1  ,rO 

—  \  v ,  .  *  —  V  .  . 

X2  1J  K 


v  *  (Cov(£  1  ))  =■  1  V 
J  X 

■  icovC  £..)) 


Thus  to  estimate  V,  comnute: 


Tq  (4): 


V.  ,  -  -  Cov(Y, ,Yj  » 


1  J  t 


X  /  ( X  —  1  ) 


Y.  Y.  -  XT.Y . 
i.w  0 ,w  i  0 


V, 
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An  example  of  a  variance  covariance  matrix  derived  in  this 
manner  is  3hown  in  Table  6-2.  Three  things  should  be  noted  when 
observing  thi3  V  matrix.  (i)  The  off  diagonal  elements  of  this 
matrix  are  definitely  not  zero.  Thi3  means  that  there  truely  is  a 
significant  covariance  relationship  as  expected.  This  will  be 
illustrated  under  (iii)  below,  (ii)  Tor  this  V  matrix,  the  main 
diagonal  elements,  which  represent  the  variance  of  the  replications 
about  the  mean  response,  are  not  significantly  different.  These 
diagonal  elements  represent  the  variances  of  the  response  about  the 
mean  of  the  various  design  points. 

The  test  used  to  determine  whether  or  not  these  variances 
are  statistically  equivalent  is  the  Burr  Foster  Q  test  statistic 
(Burr  (1974.)).  Q  is  used  here  to  distinguish  thi3  test  statistic 
from  the  9  (quadratic)  matrix. 

In  terms  of  the  sample  variances  37  computed  within  each 
treatment,  the  Burr  Foster  test  statistic  i3  given  by: 


9 


) 


2 


where  n  is  the  number  of  samples  or  design  points  in  the  optimal 
design.  Large  values  of  Q  lead  to  a  rejection  of  the  hypothesis  of 
equal  variances. 

The  calculated  Q  statistic  for  the  data  of  Table  6-2  is 
.0771 .  From  a  3urr  Foster  table  of  critical  values,  a  cf  statistic 
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of  .08  or  greater  is  considered  significant  at  the  .001  level  with 
29  degrees  of  freedom.  Thus,  the  design  points  for  the  optimal 
design  can  be  interpreted  as  having  equal  variances.  (iii)  The 
third  feature  to  note  when  examining  the  V  matrix  of  Table  6-2  is 
the  correlation  coefficient,  denoted  as  P.  This  coefficient  is 
calculated  by: 

p  m  covariance(X , Y) 
xy  (V(X)  V(Y))l/2 

It  can  be  shown  that  -1  <  P  <  1  .  The  quantity  P  is  a  measure 

of  the  association  between  the  random  variables  X  and  Y.  For 

example,  if  P  =  1 ,  X  and  Y  are  perfectly  correlated  and  the 

xy 

possible  values  of  X  and  Y  all  lie  on  a  straight  line  with  a 
positive  slope  in  the  (X,Y)  plane.  If  P  =0  the  variables  are 

xy 

3aid  to  be  unassociated,  that  i3,  linearly  unassociated  with  each 

other.  Calculating  the  value  of  P  in  this  manner  for  the  data  of 

xy 

Table  6-2,  results  in  finding  that  all  of  the  correlation 
coefficients  are  "=  .99.  It  is  quite  clear  from  this  data  that  X 
and  Y  are  nearly  perfectly  correlated  for  thi3  example.  It  is  also 
quite  clear  that  the  off  diagonal  elements  are  not  zero.  Thi3 
indicates  that  weighted  least  squares  regression  analysis  should  be 


considered . 


I.**  '•*»»*'« 
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6.'7  Weighted  Least  Squares  Regression  Analysis 

Derivation  of  the  V  matrix  w as  done  in  section  6.6  using  a 
set  up  for  weighted  least  squares  analysis.  Estimation  of  the 
benefit  function  coefficients  and  their  variance  for  weighted  least 
squares  is  as  shown  below. 

For  the  weighted  least  squares  method: 

b  =  (x,v'1xr1rv~1Y 


where  V  is  the  same  variance  covariance  matrix  as  before. 


If  the  observations  were  independent: 

2  o 


*\  o 

o  <r2 

i 


where  some  of  the  <rT  may  be  equal. 

Recall  however  that  taking  the  same  sequence  of  random 

draws  for  each  design  point  results  in  a  highly  dependent 

structure.  Consequently  the  full  variance  covariance  matrix  (V)  of 

the  observations  must  be  considered.  Thi3  V  matrix  will  be  ('5  x 

•  51  where  the  i,.jth  entry  will  be: 

Cov(Y.  , Y .  ) 
iw’  jw' 

and  the  covariance  is  defined  as: 


5q( 5 ) : 


:ov(Y.  , Y  .  )-l 
iw’  jw  k 


w  Jw  1  J 


This  is  the  same  expression  obtained  before  as  equation  4. 


An  equivalent  expression  for  this  Covariance  which  circumvents  some 
of  the  numerical  roundoff  errors  associated  with  these  calculations 


is : 

1  ( 

CT 

Given  this  V  matrix  calculated  a3  equation  (5)  or  (6),  the 
variance  of  the  coefficients  for  weighted  least  squares  analysis  is 
expressed  as : 

3q  (7): 

V(b)  -  (x-v’x)-1 

6.3  Covariance  Matrix  Singularity 

The  fact  that  there  is  a  highly  correlated  structure  due  to 
the  procedure  for  selecting  and  applying  common  random  draws  for 
all  design  points  cannot  be  ignored.  Thi3  fact  should  suggest  the 
use  of  a  weighted  least  squares  approach.  Because  of  this 
significant  correlation,  a  near  singular  V  matrix  is  derived.  For 
most  problems  this  occurs  for  the  period  T  data.  In  these  cases 
use  of  the  weighted  least  squares  formulas,  requiring  inversion  of 
the  V  matrix  simply  do  not  work. 

The  reason  for  this  singularity  is  that  due  to  the  method 
of  applying  common  random  numbers,  there  i3  a  strong  dependency 
between  responses  over  the  entire  design  3pace.  So  strong  in  fact 
that  the  correlation  coefficients  for  the  covariance  are  nearly  all 


V 


’  <F> 


155 


.99.  This  suggest  that  all  of  the  responses  lie  in  or  very  near  a 
hyperplane  in  15  space.  For  the  case  of  15  design  points,  given 
the  other  14  points,  the  15 J‘  can  be  predicted  with  very  good 
accuracy.  This  interpretation  follows  from  the  fact  that  the 
noninvertability  of  the  V  matrix  implies: 

iet(V)  •  0 

This  is  the  same  as  saying  there  exists  a  vector  "a"  such  that: 

a*7a  »  0 


which  implies: 


Var  (^>  a  Y.  )  »  0 
i«1  1  1 


This  then  comes  from  the  definition  of  a  hyperplane  in  15  space: 


Now,  a  hyperplane  in  15  3pace  can  be  thought  of  as  shown  in 
Figure  5-4.  Some  of  the  problems  start  out  as  invertable  weighted 
least  squares  problems,  however,  all  of  the  responses  fall  within 
this  hypervolume  in  15  space.  Initially,  this  hypervolume  may  be 
large  enough  or  thick  enough  that  the  7  matrix  is  invertable. 
Taking  the  result  of  the  inversion  and  estimating  the  next  set  of 
coefficients  has  the  effect  of  squeezing  down  this  hypervolume 
until  eventually  it  is  too  thin  to  allow  an  invertable  form  of  the 


V  matrix. 


Figure  6-4 

Illusion  of  Hypervolume  in  15  Spa^ 
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To  support  this  theory  a  linear  regression  wa3  performed  on 
the  response  data  using  Y1  as  the  dependent  variable  and  the 
remaining  14  mean  responses  as  the  independent  variables.  The 
model  to  be  fit  was  thus: 

The  results  were  even  more  convincing  than  anticipated.  The 
results  of  this  regression  are  shown  in  Table  5-3(a). 

Next  a  similar  regression  was  performed  on  a  set  of  data 
that  was  originally  invertable  and  observed  that  in  fact  this 
matrix  was  nearly  singular  from  the  outset.  Table  6-3(b)  shows  the 
results  of  this  regression. 

These  tables  conclusively  support  the  theory  that  the  means 
of  the  responses  all  lie  in  a  hypervolume  in  15  space.  Further, 
this  totally  explains  the  noninvertibility  of  the  variance 
covariance  matrix. 

Thus,  while  weighted  least  squares  analysis  is  indicated  in 
nearly  all  of  the  data,  its  use  brings  about  its  own  demise  due  to 
the  noninvertibility  of  the  V  matrix.  Consequently,  the 
recommended  methodology  does  not  utilize  estimation  of  the 
coefficients  using  weighted  least  squares. 

Table  6-1  in  section  6.5  showed  the  number  of  draws 
required  to  get  some  specified  accuracy  of  estimating  the  V  matrix. 
If  the  problem  did  not  change  from  period  to  period,  the  V  matrix 
could  be  estimated  one  time  using  a  large  number  of  draws,  and  then 
a  fewer  number  of  draws  could  be  used  for  subsequent  periods. 


Table  6-3 


Results  of  Linear  Regression 


(Y,  with  Y?  through  Y^) 

(a) 

(b) 

Not 

Invertible 

Invertible 

R-Squared 

1 .0 

1  .0 

Standard  Dev. 

.00835 

.1104 

Residual(SS) 

.00599 

.  1 9281 

?-Value 

354422760. 

1 400488 . 

Significance 


.0 


.000 


Since  the  same  number  of  draws  will  be  required  for  each  period 


another  criteria  for  determining  the  number  of  draws  is  to  select 
k '  so  as  to  control  the  maximum  variance  of  the  predicted  error. 
This  is  different  from  the  derivation  of  the  number  of  draws  for 
the  V  matrix  in  the  following  way.  Recall  that  in  deriving  the  V 
matrix  equation  (4)  was  used  where: 


V. 


U 


1  Cov(Y. ,Y  .) 

IT  J 


Here  K  i3  the  number  of  draws  which  is  determined  by  the  selection 
of  an  acceptable  relative  error  as  in  Table  5-1 . 

Controlling  the  maximum  variance  of  the  predicted  error 
requires  that  the  variance  covariance  terms  derived  in  equation  (4) 
be  used  but  without  the  division  by  K.  Let: 

Sq  (9): 


V  =  XV 
0 

This  is  the  desired  V  matrix  for  making  this  determination.  Vow 
for  ordinary  least  squares  it  i3  desired  to  select  k'  large  enough 
to  control  the  maximum  variance  about  the  mean  of  the  worst  design 
point,  or: 


3q  (10): 


Sax  (x (X'X)-1  (X'V  XKX'X)-1  )x')1  /2  (t  ,  df ) /JF  <  Del 

_  -  »  o  o  o  a/£^  — 


x  €.  X 
.  o 
where 


xQ  »  row  in  X  matrix  with  largest  variance 


k  *  number  of  draws 
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t  *  student  t  value  for  desired  accuracy  with  given  d.f. 

Del  ■  absolute  tolerable  error  (*  of  Y  ) 

o 

Rearranging: 

2  max  (x  (x'x)"1  (x'vxHx'xrV) 

„  _  0  0 

x0£* 

The  results  obtained  using  the  normal  least  squares 
calculations  as  applied  to  the  data  of  Table  6-2  follow.  In  this 
case,  the  initial  7  matrix  was  obtained  using  30  draws.  3y 
applying  equation  (9)  to  this  matrix  and  using  the  resul  t  (V  )  in 
equation  (11),  '  was  determined  to  be  approximately  20. 

Table  5-4  shows  other  results  given  various  settings  of 
alpha  and  the  tolerance  level  in  absolute  error. 

Table  6-4 


Number  of  Draws  for  Controlling  Maximum  Variance 


Number 
draws  for  V 

Tolerance 

level 

alpha 

setting 

Calculated 
#  draws 

30 

.05 

.01 

20 

30 

.01 

.04 

30 

100 

.01 

.05 

33 

100 

.01 

.01 

32 

?rora  this  table  it  would  seem  reasonable  to  use  either  30 
draws  (alpha  *  .04)  or  approximately  90  draws  (alpha  »  .01) 


depending  upon  the  desired  accuracy. 


These  two  choices  3eem 


reasonable  since  in  both  cases  the  number  of  draws  used  to  estimate 
V  and  the  calculated  number  of  draws  are  nearly  equivalent.  As  was 
done  in  section  6.3,  30  draws  will  be  used.  Again,  the  user  can 
use  more  draws  if  greater  accuracy  is  iesired. 

6.9  Selected  Methodology 

Given  the  statistical  information  just  presented,  it  is 
iesired  to  select  a  prudent  and  statistically  sound  heuristic  for 
the  general  problem.  Both  the  accuracy  of  the  model  and  its 
computational  efficiency  must  be  considered.  Experience  has  3hown 
that  each  three  reservoir  nonlinear  network  requires  approximately 
.0115  seconds  of  CPU  time.  Due  to  the  dynamic  nature  of  the  water 
resources  problem,  it  is  necessary  that  the  V  matrix  be  updated 
every  period.  Accordingly,  the  selection  of  the  number  of  draws 
(!C)  to  estimate  the  V  matrix  and  the  number  of  draws  (k')  for  the 
periods  to  follow  should  ideally  be  the  same. 

For  these  reasons,  it  is  felt  that  a  minimum  of  30  draws 
should  be  used  for  estimating  the  V  matrix  and  for  subsequent  use 
in  calculating  the  variance  of  the  estimated  parameters.  The 
estimation  of  the  parameters  will  be  calculated  using  ordinary 
least  squares  methods,  and  if  the  variances  of  the  parameters  is 
calculated,  it  will  be  done  using  the  full  7  matrix  in  the  manner 
of  equation  4. 

While  the  main  interest  lies  in  the  estimate  of  the 
parameters  it  is  instructive  to  know  what  levels  of  variability 


162 


exists.  3y  using  70  draws  per  level  for  all  periods  a  reasonably 
good  estimate  of  the  parameters  will  be  achieved  with  a  minimum  of 
computational  time. 

The  relative  error  of  the  standard  error  of  the  estimated 
variance  to  the  true  variance  is  25%.  Given  this  25;C  criteria, 
very  high  confidence  about  this  estimate  is  realized  as  indicated 
by  the  absolute  tolerance  level  and  the  i%  alpha  limit  from  the 
t  table  (see  Table  6-4). 

In  summary,  it  is  felt  that  30  draws  is  a  good  compromise 
between  accuracy  and  computational  time.  For  30  draws  per  level 
and  15  levels  per  period  the  computational  time  will  be  roughly  7 
seconds  per  period  for  the  three  reservoir  problem.  The  number  of 
periods  to  use  is  highly  dependent  upon  the  specific  problem  being 
considered,  its  intended  use  and  the  availability  of  data.  This 
will  be  addressed  mors  fully  in  Chapter  7. 


CHAPTER  VII 


7.  Example  Applications 

This  chapter  includes  several  example  applications  to 
demonstrate  the  feasibility  of  the  model.  Section  7.1  includes  two 
network  formulations  with  some  variations  applied  to  both.  Section 
7.2  includes  a  representative  application  to  the  Guadalupe  River 
3asin  in  Texa?  by  attempting  to  evaluate  a  proposed  new  system  to 
meet  the  demands  of  the  year  2020. 

7.1  Hypothetical  Problems 

This  section  includes  both  a  three  reservoir  and  a  four 
reservoir  example  problem.  There  are  two  versions  of  the  three 
reservoir  problem,  the  differences  being  in  the  number  of  arcs  and 
the  selection  of  the  arc  parameters. 

Figures  7-1 ,  7-2  and  7-3  represent  the  three  problems  to  be 
considered.  These  will  be  referred  to  herein  as  examples  1,  2  and 
3  respectively.  Although  the  node  structure  of  Figures  7-1  and  7-2 
are  the  same,  note  the  differences  in  the  number  of  arcs  and  in  the 
arc  parameters.  The  network  of  Figure  7-2  assigns  a  benefit  to 
having  a  minimum  amount  of  water  in  the  river  reaches.  This  may  be 
necessary  due  to  hydroelectric  concerns  or  perhaps  due  to  concerns 
for  aquatic  life.  Additionally,  the  benefits  for  supplying  demand 
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are  revised  in  this  figure  which  tend  to  place  more  value  on 
meeting  the  current  needs  than  the  benefits  for  demand  of  Figure 
7-1  . 

Reference  has  been  made  throughout  this  report  to  Figure 
4-1 .  The  network  of  Figure  7-1  is  identical  to  Figure  4-1  and  many 
of  the  experimental  results  are  used  throughout  this  report. 

For  Figures  7-1  and  7-2,  all  reservoirs  have  a  capacity  of 
25  units  of  water  with  an  allowance  for  5  additional  units  as  a 
costly  high  water  condition.  The  discretizations  are  (.5,  .7,  .9). 
For  Figure  7-3,  the  capacities  of  the  reservoirs  are  (40,  20,  30, 
25)  for  nodes  (l ,5,9,13)  respectively.  For  this  example  the  same 
discretizations  of  (.5,  .7,  .9)  are  used.  In  this  chapter,  Q  (the 
quadratic  matrix)  will  use  the  subscript  t.  Thi3  will  mean  that  if 
T  *  12,  will  be  the  assumed  quadratic  benefit  function  for 
period  T,  and  the  dynamic  programming  algorithm  will  generate  Qf 
for  t*1 1 ,  10,  9,...,*  in  this  order.  Additionally,  when  is 
used,  it  will  imply  the  existance  of  the  associated  linear  terms  of 

r 

the  full  quadratic,  in  addition  to  the  quadratic  matrix  itself. 
Using  this  notation,  for  Figures  7-1  and  7-2  is  taken  to  be  the 
negative  of  the  benefit  function  used  as  the  example  in  Chapter  5. 
Q„  for  the  example  of  Figure  7-3  uses  the  identity  matrix  for  the 
quadratic  matrix.  terms  for  these  problems  are  summarized  in 
Table  7-1  which  includes  the  quadratic  terms  and  the  associated 
linear  terms.  For  both  cases,  the  terms  were  arbitrarily 


selected . 
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Table  7-1 

Assumed  Conditions  for 


fl 

B1 

Figs  7-1 ,7-2 

-59-22 

Fig  7-3 

-15.0 
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-46 . 61 
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B3 

-39.69 

-25.0 

f3 

f2 
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-30.0 

f4 
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1 .0 

f§ 
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.52 
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In  analyzing  these  problems  it  is  interesting  to  observe 
the  effects  of  two  primary  factors,  (i)  the  nature  of  the  inflow 
data  and  (ii),  the  effect  in  the  long  run  of  the  assumed  values  for 
Q  .  These  two  issues  will  be  addressed  in  the  next  two  sections. 

7.1.1  Model  Response  To  Inflow  Data 

To  address  the  first  issue  Q^,  was  arbitrarily  selected  as 
shown  in  Table  7-1  and  the  inflow  pattern  was  varied.  Three 
variations  of  the  inflows  were  considered  as  shown  in  Figure  7-4a, 
b  and  c.  Figure  7 -4a  represents  inflows  which  are  low  in  the 
present  and  monotonieally  increase  to  higher  conditions  in  the 
future.  Since  the  algorithm  woric3  backward  from  the  future  to  the 
present  it  is  expected  that  the  current  decisions  will  not  insist 
that  water  be  stored  as  a  first  priority  since  the  model  foresees 
more  supply  in  the  future. 

Figure  7-4b  represents  the  opposite  inflow  situation  to 
that  just  discussed.  Here  there  is  more  water  in  the  present  but 
less  is  anticipated  in  the  months  to  come.  It  is  expected  in  this 
case  that  the  model  will  attempt  to  save  some  water  for  future  use. 

Finally,  Figure  7-4c  depicts  the  water  fluctuating 
throughout  the  time  horizon  going  from  wet  to  dry  then  wet  to  dry 
again.  This  may  be  more  representative  of  the  inflow  profile  over 
a  longer  time  horizon.  For  the  examples,  T  was  selected  to  be  12 
with  each  period  being  equivalent  to  one  month. 

These  inflow  characteristics  and  the  network  of  Figure  7-1 


will  be  referred  to  as  example  1-1,  1-2  and  1-3.  Similarily, 
applying  these  inflow  profiles  to  examples  two  and  three  will  yield 
three  cases  each. 

For  all  of  these  problems,  inflows  were  assumed  to  be 
distributed  normally  with  standard  deviations  running  from  50%  to 
50%  of  the  mean  for  all  reservoirs.  Actual  inflow  data  for  the 
profiles  of  Figure  7-4  are  included  in  Part  I  of  the  Appendix. 

The  negative  of  the  derived  benefit  functions  for  for 
each  of  the  example  1  cases  are  shown  in  Table  7-2.  Tables  7-3  and 
7-4  show  similar  results  for  examples  2  and  3  respectively. 

Some  comments  regarding  the  data  of  Tables  7-2, 3, 4  follow. 
There  are  two  primary  things  to  note  wien  evaluating  these  tables. 

The  first  item  of  interest  for  these  tables  deals  with  the 
magnitude  of  the  linear  terms  of  the  benefit  function  when  compared 
to  the  other  negative  costs  of  their  respective  networks.  As  an 
example,  consider  Table  7-4  and  Figure  7-3.  Consider  the  routing 
of  flows  as  the  flows  into  the  system  at  the  reservoir  nodes  are 
increased  from  an  initial  zero  level.  The  first  unit  of  flow  will 
be  routed  based  solely  on  the  linear  cost  contribution,  since  for 
the  nonlinear  arcs,  there  is  no  quadratic  contribution  to  the 
marginal  cost  at  zero  flow.  Consequently,  for  example  problem  3-1 
which  foresees  a  wet  future,  the  initial  allocation  of  water  is  to 
meet  current  demands.  This  is  seen  by  comparing  the  negative 
benefit  for  supplying  demand  with  the  negative  linear  terms  in  the 
benefit  function.  Users  6  and  14  will  first  be  supplied  water 
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until  their  minimum  demands  are  met  (saturation  of  one  of  the 
demand  arcs).  User  14  will  also  supply  its  secondary  demand  before 
any  flow  is  allocated  to  the  reservoir  for  future  storage.  Users  2 
and  6  will  compete  for  water  almost  immediately  with  their 
reservoirs. 

For  the  inflow  profile  of  Figure  7-4b,  the  future  is 
expected  to  be  dry.  In  this  case,  note  the  extreme  change  in  the 
linear  terms  of  the  nonlinear  arcs  which  indicate  a  strong  desire 
to  store  the  initial  allocation  of  water,  (likewise  for  the  inflow 
profile  of  Figure  7-4c).  This  pattern  also  holds  for  the  data  of 
Table  7-2.  Table  7-3  is  unique  in  that  the  linear  cost  terms 
assigned  to  the  nonlinear  arcs  do  not  change  much  at  all  over  the 
three  inflow  profiles.  This  is  a  result  of  placing  large  negative 
costs  on  the  river  arcs  of  Figure  7-2.  For  this  example,  the 
current  needs  dominate  the  entire  process  and  in  the  end,  all 
demands  (both  users  and  storage)  reflect  near  equal  priorities. 

The  second  item  of  interest  is  the  amount  of  quadratic 
effect  realised.  For  all  three  example  problems,  as  the  future 
supply  of  water  decreases,  the  quadratic  terms  increase. 
Similarily,  for  inflow  profile  1,  these  terms  are  very  small 
indicating  a  near  linear  situation.  Thus,  as  the  future  supply  of 
water  decreases,  interactive  forces  tend  to  arise  which  supports 
the  hypothesis  of  reservoir  interaction. 

One  final  comment  on  these  tables.  The  data  for  inflow 
profiles  1  and  2  cannot  be  compared  directly  with  profile  3  because 


V 


176 


there  is  more  water  overall  in  profile  three.  For  profile  3, 
period  1  has  a  very  large  amount  of  inflow,  as  do  the  next  few 
periods.  Thus,  the  results  for  this  profile  are  altered  by  current 
and  near  term  high  inflow  conditions.  This  could  be  what  is 
causing  the  biggest  difference  in  the  coefficients. 

As  an  additional  display  of  the  model  results,  Table  7-5 
shows  one  of  the  negative  benefit  functions  along  with  the 
calculated  standard  deviations  for  the  estimated  coefficients. 
These  standard  deviations  were  calculated  using  the  methods 
described  in  Chapter  6.  These  results  are  very  similar  for  all 
example  problems  and  indicate  that  in  fact,  30  draws  does  yield 
very  good  estimates  of  the  coefficients.  An  additional  test  was 
performed  using  this  data.  Rather  than  estimate  the  coefficients 
and  calculating  their  variance  using  the  normal  equations,  the  data 
was  provided  as  input  to  a  least  squares  regression  package 
(FIXREG).  The  disadvantage  of  sing  FIXREG  regularily  is  that  it 
requires  considerably  more  time  due  to  thsnumerous  statistics  it 
generates  and  to  its  built  in  plotting  capability.  However,  it  was 
used  on  occasion  to  assure  that  meaningful  results  were  being 
generated.  An  example  of  some  of  the  results  are  also  shown  in 
Table  7-5.  These  t  and  F  levels  of  significance  and  the  overall 
R^  term  conclusively  support  the  quadratic  as  a  valid  function  for 
the  regression.  This  example  is  highly  representative  of  the 


results  for  all  example  problems. 
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Table  7-5 

Statistical  Results  for  Coefficients 


Std.Dev 

30 

-250.9 

23.3 

B1 

-25-9 

.33 

B2 

-16.0 

.19 

B5 

-15.3 

.19 

B4 

.274 

.012 

85 

.099 

.002 

B6 

.103 

.002 

37 

.127 

.007 

38 

.114 

.008 

39 

.134 

.004 

R2  = 

.999- 

'-STAT 

"  3.7  “ 

53.2 

34.0 

35.0 

25.0 

49.5 
51  .5 
13.1 

14.5 
46.0 


F  Significance  *  5.9E*5 
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7.1.2  Operational  Use  of  Benefit  Functions 


Once  the  benefit  functions  have  been  derived  by  the  dynamic 
programming  procedure,  they  should  be  useful  in  an  operational 
context  to  make  decisions  on  the  optimal  distribution  of  available 
surface  water  in  any  given  month.  The  function  can  easily  be  saved 
by  storing  the  coefficients  of  the  quadratic  form  and  the  linear 
coefficients.  This  use  of  the  benefit  functions  is  illustrated  in 
this  section  using  the  benefit  functions  derived  for  period  1  for 
the  example  cases. 

To  use  the  benefit  functions  the  decision  maker  mu3t 
observe  his  reservoir  levels  at  the  beginning  of  the  month, 
calculate  or  observe  the  expected  inflows  and  run  the  single  period 
optimization  one  time.  The  resulting  flows  will  determine  the 
policy  for  period  1  . 

To  implement  this  activity  with  the  computer  codes  used  for 
the  dynamic  programming  procedure  a  few  changes  need  to  be  made  to 
the  data  set.  The  required  changes  are  very  simple  and  do  two 
things.  First,  they  force  the  logic  to  skip  certain  calculations 
that  pertain  to  multiperiod  and  mult-draw  problems.  Second, 
specific  changes  to  the  input  data  set  allow  the  user  to  specify 
his  desired  conditions. 

Using  the  <?1  benefit  functions  of  Tables  1-2,  3  and  4, 
several  period  1  conditions  were  assumed.  These  assumed  conditions 
represent  the  sum  of  the  observed  reservoir  levels  and  the  expected 
inflows  for  the  period  ( ie -  total  water  available  for  period  '). 
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With  these  total  inflows  the  optimum  flows  were  obtained  by  the 
network  algorithm.  These  flows  describe  the  optimum  decision  set 
for  each  condition.  Some  of  these  are  shown  in  Table  7-6. 

Specifically,  the  data  of  Table  7-6  represents  two  examples 
to  be  discussed.  The  first  two  columns  of  Table  7-6  reflect  two 
inflow  level  sets  for  period  1  as  applied  to  example  problem  1-2. 
The  third  and  fourth  columns  reflect  similar  inflow  level  sets  for 
example  1-3.  There  are  other,  perhaps  more  interesting,  examples 
which  involve  more  releases  and  transfers  of  water  between 
reservoirs,  but  these  were  selected  to  demonstrate  the  network 
optimization  technique.  In  all  these  cases  (except  the  last 
column)  all  the  water  available  at  each  reservoir  was  used  to 
supply  demand  at  the  reservoir  or  saved  in  the  reservoir.  No  water 
was  transferred  or  released. 

This  first  discussion  pertains  to  example  1-2  of  Table  7-6. 
Consider  the  first  column  with  inflows  to  the  three  reservoirs  of 
13/13/13.  Referring  to  Table  7-2  and  Figure  7-1,  it  is  clear  that 
the  linear  cost  coefficients  for  the  benefit  function  are  more 
negative  than  the  demand  cost  coefficients.  For  reservoir  1  ,  this 


means  that  the  initial  allocation  of  water  will  go  to  the 
reservoir.  The  question  is  how  much  water  will  be  stored  before 
any  demands  are  met?  As  water  is  stored  in  all  reservoirs,  the 


marginal  cost  on  the  nonlinear  arcs  for  storing  water  increases. 
At  some  point,  the  marginal  cost  will,  be  greater  than  the  marginal 
cost  for  supplying  demand.  For  reservoir  1  this  occurs  at  11.35 
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units  of  water.  At  this  point,  the  marginal  cost  for  storing  water 
in  reservoir  1  is  increased  to  -17.  This  is  the  breakpoint  for 
supplying  demand,  and  the  next  1.65  units  are  allocated  to  the 
demand  at  reservoir  1.  If  these  1.65  units  were  stored  rather  than 
used,  the  marginal  cost  would  be  greater  than  -17  and  hence,  the 
network  flows  would  not  be  optimal. 

Vith  regards  to  reservoir  2  and  demander  2,  all  13 
available  units  are  stored.  At  this  level,  the  marginal  cost  for 
reservoir  2  is  -10.24.  The  marginal  cost  must  increase  to  -10.0 
before  flow  will  be  allocated  to  demand.  Finally,  for  reservoir 
3,  the  initial  allocation  of  water  went  to  the  reservoir  since 
-15.7  <  -15.0.  However,  very  quickly,  the  marginal  cost  for  thi3 
reservoir  drops  to  -15.0  and  the  next  four  units  of  water  go  to  the 
user  (node  10).  Once  user  3  received  4  units,  it  is  again 
profitable  to  store  water  and  the  remaining  water  i3  in  fact 
stored.  Final  marginal  cost  for  reservoir  3  is  -10.1 9,  far  better 
than  supplying  an  additional  3  units  to  user  3  at  a  cost  of  -6.0 
per  unit. 

For  the  22/22/22  column  of  example  1-2,  similar  flows  will 
occur  for  the  first  13  units.  However,  with  an  additional  9  units 
of  water  available  at  each  reservoir,  user  1  receives  all  7  units 
requested  with  the  remaining  2.65  units  being  stored.  In  this  case 
the  ending  marginal  cost  for  reservoir  1  is  -13.22  versus  -17.0 


5 

4 
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from  before.  User  2  now  receives  all  3  units  demanded  at  a  cost  of 
-10.0  and  an  additional  .11  units  at  -7.0.  Given  additional  flow, 
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the  second  increment  of  demand  at  user  2  will  be  satisfied  before 
any  additional  storage.  At  the  current  levels,  the  marginal  cost 
for  reservoir  2  is  -7.0.  Reservoir  3  ending  marginal  cost  is  less 
than  -6.0  and  user  3  does  not  receive  its  second  increment  of 
demand.  Note,  in  all  these  cases,  no  flow  is  released  downstream. 
This  would  not  occur  unless  one  of  two  things  were  to  happen. 
First,  if  an  upstream  reservoir  has  an  over  abundance  of  water  such 
that  some  marginal  cost  downstream  (user  or  storage)  was 
profitable,  or  2,  if  a  reservoir  has  enough  water  to  cause  its 
marginal  cost  to  go  to  2ero  while  meeting  all  demands.  Then  it 
would  be  profitable  to  release  water  at  a  cost  of  zero  rather  than 
store  it  at  some  positive  cost  (since  additional  storage  would 
cause  the  marginal  cost  to  go  positive) . 

Turning  now  to  the  data  for  example  1-3,  for  inflows  of 
13/13/13  reservoir  1  marginal  cost  goes  to  -17.0  with  only  6. 31 
units  of  flow  with  the  next  6.69  units  going  to  user  1  at  a  cost  of 
-17.0.  This  breakpoint  is  far  lower  than  for  example  1-2  due  to 
the  significantly  higher  quadratic  terms  of  Table  7-2.  Reservoir  2 
ending  marginal  cost  is  -10.3,  still  less  than  -10.0  for  its  user. 
And  for  reservoir  3,  its  ending  marginal  cost  was  -3.4-5  (nearing 
zero) . 

These  seemingly  low  ending  reservoir  levels  are  apparently 
due  to  the  high  inflows  in  the  first  few  periods  of  profile  3. 
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For  the  18/13/18  inflows  of  example  1-3,  the  important 
thing  to  note  is  that  ending  marginal  cost  for  reservoir  3  is  zero. 
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In  this  case  rather  than  store  additional  water  (since  all  demands 
are  met),  5.972  units  are  releasd  at  zero  cost. 

Many  variations  of  these  problems  were  run  (not  always 
having  the  same  inflows  for  all  reservoirs)  with  different  results 
for  each  case. 

Some  interesting  alternatives  to  these  examples  might  be  to 
allow  transfer  of  water  back  upstream  at  zero  or  some  very  low  cost 
thereby  making  it  profitable  to  spend  a  few  dollars  to  get  the 
water  where  it  is  needed  the  most  rather  than  release  it  at  zero 
cost.  Another  interesting  alternative  would  be  to  supply  a  large 
amount  of  water  to  reservoir  t  with  zero  inflows  at  reservoirs  2 
and  3-  If  the  zero  cost  arcs  of  the  system  were  given  a  large 
capacity,  enough  water  could  be  put  into  the  system  to  meet  all 
demands  (maximize  the  current  return)  and  to  drive  the  marginal 
cost  for  all  reservoirs  to  zero  (maximize  the  future  return).  This 
could  be  done  using  slack  inflow  at  reservoir  1  at  zero  cost  and  by 
putting  a  high  penalty  on  releasing  water  to  the  ocean.  In  this 
case  the  objective  function  would  achieve  its  absolute  minimum, 
thus  deriving  the  solution  to  the  quadratic  problem. 

As  indicated  throughout  this  report,  the  reservoir  contents 
at  the  beginning  of  the  period  do  not  have  to  correspond  to  the 
discretized  reservoir  water  levels  used  to  derive  the  benefit 
functions.  The  functional  form  of  the  benefit  function  was  derived 
from  the  discretized  levels,  but  this  form  is  a  continuous  form  and 
applies  for  any  set  of  water  level  combinations  that  fall  within 
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t  7.1.3  Effect  of  on  the  Model 

The  main  concern  here  ie  the  effect  that  the  assumed  values 
for  Q,j  have  on  the  model.  Naturally,  all  other  arc  cost  parameters 
are  also  questionable  since  selecting  them  is  an  art  in  itself. 

*  However,  it  is  assumed  that  the  user  will  have  a  fairly  good  feel 

for  these  values.  Vhat  he  will  not  have  a  good  feel  for  is  the 
future  value  of  water  at  the  end  of  the  time  horizon. 

To  measure  this  effect,  the  network  and  inflows  for  example  \ 

1-2  were  used.  Several  variations  of  the  period  T  quadratic 

benefit  functions  and  linear  terms  were  examined  with  nearly 
identical  results  in  all  cases.  The  results  3hown  in  Table  7-7 
represent  three  of  these  conditions.  This  table  shows  all  12  of 
the  Qt  benefit  functions  as  they  were  derived  over  time.  For  the 
three  cases  shown,  case  1  started  with  the  assumed  benefit  function 
of  Chapter  5.  Case  2  kept  the  same  linear  terms,  but  used  the 

* 

identity  matrix  as  the  quadratic  matrix.  Case  3  used  the  same  S 

quadratic  matrix  as  case  1  but  changed  the  linear  terms.  Observing 
the  first  few  periods  of  results  (periods  12, 11, etc.),  it  is  quite 

* 

clear  that  they  are  significantly  different.  However,  it  is  noted 
»  that  after  approximately  8  periods  of  data  (Q^)  the  coefficients  of 

the  benefit  functions  are  very  close  in  value.  After  12  periods, 
all  three  Q1  benefit  functions  are  nearly  identical.  Many  other 
variations  were  tried  with  similar  results.  Consequently,  it  is 
determined  that  for  these  problems,  a  minimum  of  3-12  periods  are 
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Table  7-7 


« 


Period 


3 


1 1 


0 


10 


Benefit  Function  Convergence 


Coeff. 

Case  1 

Case  2 

J.aas.J_ 

BO 

0.0 

0.0 

0.0 

B1 

-59.22 

-59-22 

-50.22 

B2 

-46.61 

-46.61 

-40.61 

B3 

-39-69 

-39.69 

-30.69 

B4 

.3633 

1  .0 

1 .0 

B5 

.5320 

1  .0 

1  .0 

B6 

.5226 

1  .0 

1  .0 

B7 

.6340 

0.0 

0.0 

BS 

.4066 

0.0 

0.0 

39 

.6770 

0.0 

0.0 

BO 

-571 . 

-351  . 

-475. 

B1 

-25.38 

-49-11 

-22.68 

B2 

-25.41 

-32.10 

-21 .37 

B3 

-24.99 

-31.25 

-21 .67 

B4 

.155 

.770 

.147 

B5 

.248 

.430 

.317 

36 

.313 

.490 

.437 

37 

.111 

.01 1 

.000 

B3 

.076 

.000 

.000 

39 

.287 

.263 

.081 

30 

-203. 

-316. 

-179. 

B1 

-19.72 

-34.44 

-13.93 

32 

-19.94 

-27.00 

-17.53 

B3 

-22.46 

-23.18 

-17.53 

B4 

.048 

.420 

.046 

B5 

.154 

.270 

.181 

B6 

.266 

.350 

.271 

B7 

.031 

.004 

.000 

B8 

.026 

.001 

.000 

39 

.208 

.291 

.098 

t 
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Table  7-7  (continued) 
page  2 


30 

-171  . 

-275. 

-153. 

31 

-13.26 

-25-56 

-17.90 

B2 

-17.89 

-23-71 

-15.71 

B3 

-20.15 

-24.37 

-17.93 

B4 

.023 

.21 1 

.023 

B5 

.116 

.190 

.121 

B6 

.201 

.253 

.214 

37 

.012 

.002 

.000 

B8 

.01 1 

.001 

.000 

B9 

.184 

.237 

.108 

30 

-134. 

-172. 

-1  23. 

B1 

-17.68 

-21  .15 

-17.67 

32 

-16.77 

-22.52 

-14.71 

33 

-19.14 

-23.20 

-16.68 

B4 

.014 

.101 

.017 

35 

.091 

.161 

.083 

B6 

.139 

.199 

.156 

B7 

.005 

.001 

.000 

38 

.005 

.001 

.000 

B9 

.165 

.290 

.112 

BO 

-208. 

-260. 

-197. 

B1 

-17.52 

-19.58 

-17.25 

B2 

-15.47 

-20.06 

-14.01 

33 

-16.36 

-20.63 

-15.35 

B4 

.013 

.041 

.007 

35 

.079 

.133 

.073 

B6 

.113 

.157 

.128 

B7 

.003 

.001 

.001 

38 

.002 

.001 

.000 

39 

.141 

.241 

.104 

4 

i 
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Table  7-7  (continued) 
page  3 


BO 

-207. 

-237. 

-i  97. 

81 

-21 .01 

-21  .19 

-20.94 

32 

-14.96 

-17.49 

-14.24 

33 

-15.29 

-17.75 

-14.76 

34 

.114 

.119 

.  114 

35 

.077 

.106 

B6 

.096 

.120 

.105 

37 

.029 

.028  ■ 

.027 

B8 

.022 

.023 

.020 

B9 

.136 

.  1  98 

.119 

BO 

-204. 

-226. 

-199. 

B1 

-24.48 

-24.59 

-24.39 

B2 

-15.34 

-16.37 

-14.34 

33 

-15.53 

-17.06 

-15.16 

B4 

.219 

.220 

.220 

B5 

.081 

.098 

.082 

B6 

.097 

.111 

.106 

37 

.063 

.064 

.061 

B8 

.054 

.056 

.052 

B9 

.147 

.186 

.136 

BO 

-223. 

-236. 

-220. 

B1 

-23.61 

-23.73 

-23.52 

B2 

-14.65 

-15.58 

-14.32 

B3 

-14.38 

-15.79 

-14.67 

34 

.208 

.211 

.210 

B5 

.075 

.084 

.076 

B6 

.088 

.095 

.096 

37 

.058 

.057 

.055 

B8 

.052 

.053 

.049 

B9 

.140 

.163 

.134 
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Table  7-7  (continued) 
page  ^ 


30 

-1^6. 

-1  36. 

-174. 

31 

-25.38 

-26.0 

-25.36 

32 

-16.67 

-17.14. 

-16.43 

33 

-16.52 

-16.96 

-16.34 

B4 

.248 

.252 

.252 

35 

.101 

.105 

.103 

36 

.111 

.113 

.117 

37 

.124 

.124 

.  1  22 

38 

.096 

.097 

.092 

39 

.179 

.190 

.175 

30 

-215. 

-221  . 

-213. 

B1 

-23.30 

-29.33 

-29.31 

32 

-13.79 

-19.04 

-19.59 

B3 

-19.03 

-19.29 

-17.95 

34 

.322 

.322 

.331 

35 

.140 

.142 

.141 

36 

.141 

.142 

.149 

37 

.195 

.195 

.193 

38 

.  1 46 

.146 

.144 

39 

.232 

.238 

.235 

30 

-249. 

-250. 

-247. 

B1 

-25.92 

-25-92 

-25.92 

32 

-15.94 

-16.00 

-15.  -73 

33 

-1 5.74 

-15-92 

-15.57 

34 

.274 

.274 

.273 

35 

.098 

.099 

.098 

B6 

.102 

.103 

.104 

37 

.129 

.127 

.126 

38 

.114 

.114 

.111 

39 

.191 

.184 

.178 
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required  to  dampen  out  the  effect  of  the  assumed  benefit  function 
for  period  T 

v 

7.1.4  Computation  Times 

This  section  includes  a  brief  summary  of  computation  times. 
Table  7-8  shows  the  CPU  times  and  I/O  times  required  for  these 
problems.  All  problems  were  run  on  the  University  of  Texas  CDC 
6600  system.  The  results  indicate  that  the  algorithm  requires 
approximately  .0115  seconds  of  CPU  time  for  each  three  reservoir 
nonlinear  network  and  approximately  .028  seconds  for  each  four 
reservoir  network. 

7.1.5  Benefit  Contours 

3efore  going  on  to  the  application  to  the  Suadalupe  River 
Basin,  it  is  interesting  to  observe  certain  properties  of  the  model 
as  it  relates  to  the  hypothetical  examples.  Specifically,  for  the 
three  reservoir  problems,  it  is  possible  to  plot  iso-benefit  curves 
as  a  function  of  two  reservoirs  while  holding  the  third  reservoir 
at  some  specified  level.  This  may  be  important  for  some  systems 
that  for  one  reason  or  another  require  that  a  given  reservoir  be 
closely  regulated  and  maintained  at  or  near  a  precise  level.  Any 
of  the  three  reservoirs  can  be  specified  as  fixed  for  the 
algo rithm. 

As  an  example,  suppose  it  is  of  interest  to  fix  reservoir 
three  at  20  units  of  water  for  example  1-2.  In  so  doing  it  i3 
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Table  7-8 


Summary  of  CPU 

and  I/O  Times  (s 

econds) 

Problem 

CPU  Time  I/O  Time 

1  -1 

58.3 

10.6 

1-2 

70.6 

11  .0 

1-3 

64.9 

10.5 

2-1 

63.3 

10.5 

2-2 

54.0 

10.7 

2-3 

60.2 

10.6 

3-1 

258.9 

12.6 

3-2 

267.5 

12.6 

3-3 

246.4 

12.5 

Network  Solution  Times 

Number  of  Number 

Reservoirs  Networks 

of  Average 

Solved  CPU  Time 

Time  Per 
Network 

3 

5400  62.23 

.0115 

4 


9000 


257,6 


.023 


desired  to  observe  the  resulting  iso-benefit  curves. 

A  program  (CONTI )  has  been  written  which  takes  the  desired 
benefit  function  and  target  total  benefits  as  input  data.  This 
program  then  calculates  and  plots  the  desired  curves  of  equal 
benefit.  Figure  7-5  shows  a  typical  plot  for  this  problem.  In 
this  figure  it  is  observed  that  a  major  part  of  the  ellipses  lie 
outside  the  feasible  region  determined  by  the  reservoir  capacities. 
There  are  an  infinite  number  of  these  ellipses,  depending  upon  the 
target  total  benefit.  For  Figure  7-5,  the  center  of  the  ellipse  is 
far  outside  the  feasible  region.  Recall  that  for  this  problem,  all 
reservoirs  were  restricted  to  a  maximum  capacity  of  25  unit3  with  5 
additional  units  allowed  at  a  penalty.  Thus,  these  curves  ar9  only 
meaningful  in  the  region  of  30  units  or  less  for  each  reservoir. 
If  the  feasible  region  were  unbounded,  the  center  of  the  ellipse 
would  represent  the  maximum  future  benefit  given  reservoir  3  (f^) 
is  held  at  20  units.  If  f^  were  not  fixed,  the  true  optimal 
solution  for  this  benefit  function  could  be  obtained  using  the 
fo llowing: 

-  0  for  i  -  1  ,2,3 
i 

The  results  for  this  problem  are  f1  *  34.5,  f2  ■  29.6  and  f^  •  32 
with  a  maximum  benefit  of  -1140. 

Now  letting  f^  *  32,  a  similar  plot  of  benefit  functions  is 
obtained  as  shown  in  Figure  7-6.  This  plot,  as  expected,  has  its 
center  at  (34.5,  29.6)  with  the  desired  benefit  of  -1140.  This 

point  indicates  the  maximum  of  the  benefit  function.  The  user 
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Figure  7-6  Optimum  Ill ipse.  Future  Only 
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would  not  operate  here  however  since  he  must  also  consider  the 
tradeoffs,  with  the  current  returns.  Also,  the  function  is 
certainly  not  valid  outside  of  the  feasible  region. 

As  mentioned  before,  some  of  these  ellipses  fall  mostly 
outside  the  feasible  region.  It  is  somewhat  more  useful  to 
concentrate  on  the  feasible  region  only.  Figure  7-7  has  its  range 
limited  to  a  maximum  value  of  30.  For  the  data  of  example  1-2, 
only  parts  of  the  ellipses  are  contained  in  this  region,  and  the 
maximum  benefit  attained  is  at  the  upper  right  hand  corner  of  the 
feasible  region.  In  this  case,  f  1  =  30,  f^  ”  30  and  f^  =*  20.  The 
value  of  the  benefit  function  at  this  point  is  -1112.  This  is  as 
expected  since  this  example  attempts  to  force  the  storage  of  water 
for  the  future.  Given  a  fixed  amount  of  water  in  the  system,  it  is 
unlikely  that  these  levels  would  be  stored  due  to  the  tradeoff  of 
current  requirements  and  the  fact  that  to  reach  levels  of  30  in  any 
of  the  reservoirs  would  require  that  penalty  arcs  be  used.  While 
this  is  allowed  and  could  be  profitable,  it  would  depend  on  the 
costs  assigned  to  other  network  arcs. 

Another  set  of  plots  taken  from  a  different  trial  problem 
(and  hence  a  different  benefit  function)  had  the  form  of  Figure 
7-8.  The  center  of  the  ellipses  are  contained  in  the  feasible 
region.  This  type  of  curve  will  occur  when  the  model  determines 
that  it  is  not  too  important  to  save  water  for  the  future,  (example 
1-1  data). 

One  interesting  thing  to  note  regarding  these  ellipses  is 
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the  amount  of  rotation  that  they  have.  This  rotation  is  an 
indication  that  some  interaction  between  the  plotted  reservoirs  is 
present.  If  no  interaction  were  present,  the  main  axes  of  the 
ellipses  would  parallel  the  x  and  y  axes. 
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7.2  Guadalupe  River  Basin  -  Stochastic  Case 

The  Guadalupe  and  San  Antonio  River  Basins  are  situated  in 
the  southern  Coastal  Bend  area  of  Texas.  From  their  headwaters  in 
the  Edwards  Plateau  region  of  Central  Texas,  the  two  rivers  cross 
the  Gulf  Coastal  Plains  and  combine  into  a  single  stream  shortly 
before  entering  San  Antonio  Bay.  The  average  annual  rainfall  over 
the  two  drainage  basins  varies  from  35  inches  near  the  Gulf  Coast 
to  25  inches  in  the  area  of  the  headwaters. 

The  water  needs  in  the  basins  are  presently  supplied 
largely  from  groundwater  formations  underlying  the  region,  with  the 
principle  source  of  groundwater  being  the  Edwards  and  associated 
Limestone  Aquifer.  This  underground  reservoir  i3  a  water  supply 
source  for  irrigation  and  for  many  municipalities,  including  the 
city  of  San  Antonio.  The  Edwards  Aquifer  is  also  the  source  of 
water  for  the  Comal  and  San  Marcos  Springs.  These  springs  provide 
the  major  portion  of  base  flow  in  the  Guadalupe  River. 

Canyon  Reservoir  is  the  only  existing  major  storage 
reservoir  in  the  Guadalupe  River  3asin.  The  reservoir  provides 
both  water  supply  and  flood  control  storage.  Six  small 
hydro-electric  dams  on  the  Guadalupe  River  downstream  from  Sew 
Braunfels  constitute  the  only  other  significant  impoundments  in  the 
Guadalupe  3asin. 

The  consumptive  water  requirements  within  the  Guadalupe 
3asin  are  currently  being  adequately  supplied  from  a  combination  of 


surface  and  subsurface  sources.  However,  with  the  greatly 
increased  water  demands  projected  for  this  area,  it  is  evident  that 
judicious  water  resources  management  will  be  essential  in  order  to 
fully  supply  future  water  needs. 

The  major  demand  center  for  these  two  basins  is  the  San 
Antonio  metropolitian  area.  The  city  and  adjoining  suburbs  have 
experienced  rapid  population  growth  in  the  past  several  decades, 
and  are  projected  to  have  greatly  increased  populations  in  the 
future.  The  area’s  water  needs  in  the  past  have  been  supplied 
exclusively  by  pumpage  from  the  Edwards  Aquifer,  however,  this 
municipal  and  industrial  pumpage  added  to  the  irrigation  pumpage  in 
the  Balcones  Escarpment  area  to  the  west  of  the  city  is  presently 
approaching  the  average  annual  recharge  into  the  Edwards  formation. 
Should  the  increased  future  irrigation,  municipal  and  industrial 
water  requirements  continue  to  be  supplied  from  groundwater,  the 
Edwards  Aquifer  will  undergo  drastic  reductions  in  water  levels, 
diminishing  springflows  and  severe  deterioration  in  water  quality. 

As  an  alternative  to  this  depletion  of  the  Edwards  Aquifer, 
the  Texas  Water  Development  Board  (1975)  has  proposed  the 
development  of  a  conjunctive  ground  and  surface  water  resources 
system  to  supply  the  San  Antonio  area.  The  plan  calls  for  limiting 
San  Antonio  pumpage  to  215,000  acre-feet  annually,  and  supplying 
the  remaining  water  requirements  with  surface  water  conveyed 
through  pipelines  from  a  system  of  reservoirs  in  the  San  Antonio 
and  Ouadalupe  River  Basins.  This  application  will  evaluate  the 


proposed  changes  to  the  Guadalupe  River  Basin. 


A  number  of  reservoir  and  pipeline  projects  were 
considered  as  possible  components  in  the  optimal  water  supply 
system.  Of  concern  here  is  the  addition  of  three  reservoirs  to  go 
along  with  the  existing  Canyon  Reservoir.  These  three  new 
reservoirs  include  Cloptin  Crossing,  Cuero  I  and  Cuero  II.  A 
schematic  of  the  proposed  3ytem  was  shown  in  Figure  3-16. 

The  deterministic  solution  to  the  Guadalupe  River  Basin 
application  was  discussed  in  Chapter  3  and  shown  in  the  Appendix. 
The  deterministic  inflows  were  taken  to  be  the  mean  values  of  the 
inflows  by  month  over  the  4-6  year  period.  The  demands  represented 
the  projected  incremental  demands  over  the  1970  requirements. 
Total  demands  are  not  used  since  this  system  was  designed  to 
satisfy  the  demand  growth,  with  the  assumption  that  the  current 
methods  for  supplying  demand  would  continue  to  be  used  at  their 
current  levels.  In  this  case  it  was  shown  that  sufficient  water 
would  be  available  on  the  average  to  meet  the  projected  demands  for 
the  year  2020.  3ased  upon  this  result,  it  is  expected  that  the 
stochastic  solution  might  also  indicate  that  sufficient  water  will 
be  available. 

The  stochastic  problem  faces  many  different  situations  than 
the  deterministic  problem.  The  two  most  important  of  these  are  the 
fact  that  inflows  can  vary  considerably  from  the  deterministic  case 
and  that  due  to  the  uncertainties  involved,  the  storage  of  water  is 
a  much  more  significant  factor.  Recall  that  in  the  deterministic 
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*  case,  the  reservoirs  tended  to  be  held  at  their  minimum  levels. 
For  the  stochastic  problem,  this  would  generally  not  be  the  case. 

♦ 

,  The  raw  data  for  this  problem  was  supplied  by  the  Texas 

Water  Development  Board  (1975)  which  included  the  proposed 
reservoir  design,  the  reservoir  capacities,  the  benefits  for 
meeting  the  specified  demands  and  the  monthly  inflow  data. 

The  formulated  model  is  3hown  in  Figure  7-9-  Note  that 
this  is  exactly  the  same  model  as  each  of  the  12  periods  for  the 
deterministic  case. 

As  in  the  deterministic  case,  the  total  annual  demand  will 
be  spread  evenly  over  the  12  months  of  the  multiperiod  problem. 
Thi3  is  reflected  in  the  arc  capacities  going  from  the  reservoirs 
or  junction  nodes  to  the  demand  nodes. 

Table  A-3  in  Part  II  of  the  Appendix  lists  the  mean 
inflows,  which  were  used  for  the  deterministic  problem,  along  with 
their  calculated  standard  deviations.  Other  data  for  the  Guadalupe 
River  3asin  is  also  shown  in  Part  II  of  the  Appendix.  For  the 
stochastic  example  problems  in  section  7-1 ,  the  inflow  data  was 
assumed  to  be  distributed  normally.  However,  a  brief  survey  of  the 
data  for  this  problem  reveals  that  the  data  is  most  likely  not 
normal.  A  Shapiro-Wilk  (1965)  test  for  normality  wa3  conducted 
which  verified  this  conclusion.  The  Shapiro-  Wilk  V  statistic  was 
calculated  to  be  .872.  Thi3  value  would  have  to  be  .924  or  greater 
to  accept  the  normal  hypothesis  at  the  .01  level  of  signif icance. 


Many  water  resources  managers  are  concerned  with  the 
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occurances  of  floods  and  their  analysts  typically  consider  the 
inflow  distribution  for  floods  to  be  log-normal  or  in  some  oases 
log-normal  Pearson  type  III.  Although  this  problem  is  concerned 
with  the  operation  of  a  system  of  reservoirs  under  all  inflow 
conditions,  it  was  decided  to  test  the  data  to  see  if  it  was 
distributed  log-normally.  This  test  had  to  be  altered  some  since 
there  were  several  months  whose  inflows  were  zero  which  is  not 
allowed  for  the  log-normal  case.  Accordingly,  the  2ero  inflows 
were  ignored  and  a  test  for  normality  was  put  to  the  logarithms  of 
the  remaining  data.  Again  the  results  indicated  that  this  data  did 
not  possess  the  characteristics  of  a  log-normal  distributor  The 
ohapiro-Wilk  ¥  statistic  was  calculated  to  be  .3629. 

After  several  other  considerations  were  ruled  out,  it  was 
decided  to  simply  use  the  emperical  data  as  is,  and  supply  the 
random  number  sequence  to  corresponding  inflows.  Since  the  data 
was  already  arranged  for  all  four  reservoirs  by  month  and  year,  the 
question  of  correlating  the  inflows  between  reservoirs  was  taken 
care  of.  In  the  earlier  examples  the  3ame  random  number  was 
applied  to  all  reservoirs.  Thus,  the  inflows  for  all  reservoirs  in 
the  basin  were  perfectly  correlated.  For  this  problem,  if  the 
random  number  for  the  month  of  July  turns  out  to  be  33,  the  inflows 
for  all  reservoirs  will  be  determined  by  selecting  the 
corresponding  July  inflows  in  the  year  1962  ( 1 925  *38-1).  The 
random  number  generator  will  generate  a  number  between  1-46  from  a 
uniform  distribution.  This  number  will  represent  the  desired  row 


in  the  inflow  matrix  for  each  of  the  reservoirs. 

One  thing  to  note  with  regards  to  this  data  is  the  extreme 
fluctuation  of  inflows.  The  most  severe  fluctuation  occurs  at  the 
Cuero  I  reservoir  where  in  the  month  of  July,  the  range  of  inflows 
varied  from  zero  acre- feet  in  1963  and  1964  to  648  acre-feet  in 
1936.  Other  typical  fluctuations  range  from  a  low  of  zero  to  200+ 
acre-feet  for  some  months.  This  factor  alone  indicates  the  severe 
instability  and  'inreliability  of  water  supplies  for  this  area  of 
Texas. 

3ased  upon  the  decisions  made  in  Chapter  6,  this  model  was 
run  for  12  periods  using  30  draws  per  level.  The  discretized  range 
of  levels  was  selected  to  be  (.2,. 55,. 9).  For  this  four  reservoir 
model  there  will  be  25  level  combinations  in  the  experimental 
design  matrix.  At  30  draws  per  level,  there  will  be  750  nonlinear 
network  problems  to  solve  for  each  period  or  a  total  of  9000  for 
the  12  periods.  The  assumed  for  this  problem  were  as  shown  in 
Table  7-9  (col  a).  The  coefficients  for  the  Q1  benefit  function 
derived  using  the  dynamic  programming  algorithm  are  shown  in  Table 
7-9  (col  b) 

Before  discussing  the  results  of  this  application  a  few 
points  should  be  made.  First  of  all,  only  46  years  of  data  were 
available.  Ideally,  this  is  not  enough  data  to  attempt  to 
characterize  the  distribution  of  inflows  and  hence  supports  the  use 
of  the  emperical  data.  Secondly,  the  random  selection  of  30  draws 
was  done  with  replacement.  Finally,  observing  Table  A-4  of  the 
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and  31 

Table  7-9 

For  Ouadalupe  3aein 

(All  Years) 

a 

b 

*1 

31 

-500. 

-196. 

32 

-500. 

-388. 

33 

-1000. 

-82. 

34 

-1000. 

-86. 

35 

1 .0 

.16 

36 

1  .0 

1  .32 

37 

1  .0 

.on 

38 

1  .0 

.016 

39 

0.0 

.014 

310 

0.0 

.019 

31 1 

0.0 

.01 1 

31  2 

0.0 

.014 

313 

0.0 

.007 

314 


0.0 


.010 
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Appendix,  the  projected  demands  for  this  problem  are  extremely  low 
when  compared  to  the  quantities  of  water  available.  A  strict 
application  of  this  model  could  be  expected  to  yield  sufficient 
water  to  supply  all  demands. 

Returning  to  the  results  for  this  problem,  note  the  benefit 
function  for  9^  is  based  on  the  assumption  that  t  *  1  is  May. 
Thus,  the  last  month  of  the  time  horizon  is  April.  Since  it  takes 
3-1 1  months  for  the  effects  of  the  initial  conditions  to  subside, 
these  results  are  good  for  only  about  two  months.  However,  if  this 
model  were  run  for  24  months,  the  benefit  functions  for  each  of  the 
last  12  months  would  be  valid  representations  of  the  future  value 
of  water. 

The  results  of  Table  7-9  indioate  that  for  Canyon  and 
Clopti.n  Crossing,  the  initial  allocation  of  water  would  be  to 
storage  since  their  linear  terms  are  less  than  any  of  the  demand 
cost  functions.  ?or  the  two  Cuero  reservoirs,  the  top  priority  is 
to  meet  demand.  Only  Canyon  and  Cloptin  Crossing  have  significant 
terms  in  the  quadratic  function.  As  more  and  more  water  becomes 
available  these  terms  will  eventually  reduce  their  marginal  cost  to 
a  level  where  supplying  demand  will  be  profitable.  These  results 
are  expected  since  using  all  inflow  data  and  30  draws,  the  expected 
inflows  will  approach  the  mean  inflows  used  in  the  deterministic 
case.  The  design  of  these  reservoirs  was  such  that  they  would  be 
expected  to  meet  all  demands  over  a  10  year  low  flow  condition 
given  that  they  were  full  at  the  begining  of  this  10  year  period. 
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Thus,  the  above  results  using  expected  inflows  are  not  surprising 
at  all- 

in  order  to  investigate  a  particularly  dry  period  one  may 
impose  upon  the  system  a  10  year  low  flow  profile  for  inflows. 
This  wa3  done  using  the  years  1947-1956  as  the  10  year  low  flow 
period.  The  resulting  benefit  function,  Q1 ,  for  May  and  again  for 
August  are  shown  in  Table  7-10.  For  the  low  flow  profiles,  all  of 
the  linear  terms  of  the  quadratic  are  higher  than  they  were  when 
the  data  for  the  entire  46  years  were  used.  As  before.  Canyon  and 
Cloptin  Crossing  have  the  lowest  linear  cost  terms  with  the  two 
Cuero  reservoirs  marginal  costs  becoming  nearly  equal  to  the  demand 
marginal  costs.  Also,  the  terms  of  the  quadratic,  while  slightly 
more  pronounced,  are  similar  to  the  all  years  results.  Again,  for 
both  cases,  the  results  appear  to  indicate  that  the  design  is 
sufficient  to  satisfy  the  10  year  low  flow  profile.  This  is  not 
surprising  since  the  design  was  selected  to  meet  this  criteria. 

The  main  conclusion  to  be  drawn  from  this  application  is 
that  the  reservoirs  seem  to  be  adequately  designed  (if  not  overly 
designed)  to  satisfy  their  intended  purpose  of  safegarding  against 
shortages  in  the  year  2020.  The  derived  benefit  functions  could  be 
used  to  make  operational  decisions  if  this  system  were  built. 

7.3  Number  and  Duration  of  Time  Periods 

The  selection  of  the  number  and  length  of  time  periods  is 
highly  dependent  upon  the  specific  problem  being  considered. 
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Table  7-10 

Results  for  May  and  August  (Guadalupe  Basin) 


*1 

Way 

*1 

August 

B1 

-340. 

-353. 

B2 

-536 . 

t 

oi 

O 

B3 

-106. 

-115. 

B4 

-98. 

o 

1 

B5 

.332 

.341 

36 

i  .70 

1  .74 

37 

.017 

.019 

B8 

.017 

.019 

39 

.012 

.013 

310 

.023 

.030 

31 1 

.012 

.013 

B12 

.010 

.014 

313 

.005 

.006 

B14 

.012 

.013 

V. 
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Consequently,  this  section  is  intended  to  point  out  some  of  the 
possibilities  and  concerns  that  the  user  should  be  aware  of.  The 
sample  problems  in  this  chapter  dealt  with  some  of  the  specifics 
alluded  to  in  this  section. 

It  first  must  be  determined  what  types  of  decisions  the 
user  needs  to  make.  Is  he  concerned  with  long  term  or  short  term 
information?  In  the  long  term  case  he  may  want  to  let  a  time 
period  be  3  months  or  as  long  as  a  year  and  consider  a  time  horizon 
of  20  years.  This  might  be  the  case  if  he  were  concerned  with 
using  the  results  of  the  model  for  the  evaluation  of  the  design  of 
a  proposed  new  system. 

On  the  other  hand,  if  he  is  more  concerned  with  the  current 
operation  of  a  system,  he  may  want  to  specify  his  time  period  as  a 
month  or  perhaps  a  week  and  run  his  model  for  from  1  to  5  years. 
In  either  case,  he  may  also  be  limited  by  the  available  data. 
Specifically,  he  will  require  access  to  historical  runnoff  data 
which  may  not  be  available  for  a  given  choice  of  time  duration  or 
ho rizon. 

These  two  options  (and  there  are  several  others)  require 
different  sets  of  data  and  assumptions.  For  instance,  in  the  case 
of  current  operations,  he  would  be  much  more  concerned  with  the 
accuracy  of  the  model  parameters  such  as  the  benefits  for  supplying 
demand  and  would  require  better  knowledge  of  the  inflow  statistics. 

Another  concern  might  depend  upon  the  nature  of  droughts  in 
the  basin  of  concern.  Some  parts  of  the  country  are  more  sensitive 
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to  short  tern  droughts  due  to  their  immediate  effect  on  agriculture 
or  livestock  whereas  others  may  experience  droughts  which  may  last 
for  several  years. 

In  any  case,  he  would  like  to  select  a  time  horizon  that  is 
far  enough  into  the  future  such  that  the  initial  estimates  of  the 
quadratic  benefit  function  for  period  T  has  little  effect  upon  the 
functions  derived  for  the  early  periods.  For  most  of  the  trial 
data,  and  for  the  sample  problems  a  duration  of  one  month  was 
considered  and  the  time  horizon  was  restricted  to  one  year. 
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CHAPTER  VIII 
3.  Summary 

A  review  of  the  literature  indicates  that  although  much 
work  has  been  lone  in  determining  optimal  water  resources 
management  policies,  only  in  the  last  5-10  years  have  authors  begun 
to  address  larger  systems  to  include  stochastic  inflows  and 
interactive  benefits  realized  for  a  multi reservoir  system.  This  is 
important  for  many  reasons.  First  of  all  of  primary  concern  in  the 
development  of  a  realistic  water  resources  model  is  the  inclusion 
of  the  future  uncertainties  of  water  supplies.  Because  of  this 
uncertainty,  current  decisions  regarding  the  storage  or  release  of 
water  have  a  direct  effect  on  future  operations.  Provisions  must 
be  allowed  for  continued  operations  in  the  event  of  future 
shortages  or  oversupply. 

To  evaluate  the  system  given  future  uncertainties,  it  is 
felt  that  through  the  optimal  operation  of  a  multireservoir  system, 
greater  benefits  can  be  realized  than  through  the  individual 
optimization  of  each  reservoir.  This  joint  benefit  stems  from  the 
ability  to  transfer  water  between  some  reservoirs  in  the  system  and 
creates  the  need  for  a  means  of  representing  and  evaluating  these 
interactive  effects.  This  is  handled  by  generating  a  nonlinear 
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nonaeparable  quadratic  representation  for  the  future  value  of 
water. 

One  of  the  principle  contributions  of  this  research  has 
been  the  formulation  and  computer  implementation  of  an  optimization 
algorithm  for  solving  the  convex  nonlinear,  nonaeparable 
minimization  problem  for  the  generalized  network.  It  is  believed 
that  this  work  represents  the  first  application  of  nonlinear 
generalized  network  codes  modified  to  handle  such  problems,  without 
reverting  to  a  piecewise  linear  approximation.  This  modification 
builds  upon  the  highly  efficient  linear  codes  of  Jensen  and  Barnes 
(1980)  and  extends  the  practical  application  of  network  codes  to  a 
new  3et  of  problems.  The  only  restriction  for  this  problem  is  that 
the  overall  objective  function  which  is  a  combination  of  several 
linear  terms  and  some  quadratic  terms  be  a  convex  cost  function. 

Another  major  contribution  of  this  research  is  the 
formulation  and  computer  implementation  of  a  stochastic  dynamic 
programming  algorithm  for  solving  multiperiod  decision  problems 
with  uncertainty.  This  algorithm  combines  the  techniques  of 
dynamic  programming,  network  flow  programming  and  regression 
analysis  in  a  unique  way.  A  principle  feature  i3  the 
representation  of  the  value  of  the  recursive  function  as  a 
nonlinear  function.  This  functional  representation  greatly 
relieves  the  dimensionality  problems  usually  associated  with  large 
dynamic  programming  problems.  This  benefit  function  is  derived  by 
using  the  network  programs  to  optimize  the  model.  These  optimal 
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results  are  then  used  in  a  least  squares  regression  model  which 
fits  a  full  guadratic  to  the  data.  This  benefit  function  is  then 
translated  into  network  parameters  and  these  become  known  values 
for  the  solution  of  the  network  in  the  next  dynamic  programming 
period.  By  following  the  typical  dynamic  programming  recursion 
methods,  the  benefit  functions  for  each  period  are  ultimately 
generated  and  can  then  be  used  to  evaluate  decisions  in  an 
operational  context. 

Using  the  above  mentioned  algorithm,  these  procedures  are 
applied  to  a  water  distribution  problem.  Several  example 
applications  using  hypothetical  systems  are  evaluated  to  verify  and 
demonstrate  the  approach.  An  application  is  then  made  to  the 
Ouadalupe  River  Basin  in  Texas.  This  application  involves  a 
proposed  four  reservoir  system  for  this  basin  designed  to  satisfy 
the  projected  demands  of  the  year  2020.  The  results  showed  the 
Texas  Water  Development  Board  four  reservoir  design  to  be  more  than 
adequate  to  meet  the  projected  demands. 

One  of  the  primary  advantages  of  this  model  formulation  for 
the  water  resources  application  is  the  ease  of  providing  the 
required  data.  Since  the  multiperiod  model  is  simply  multiple 
copies  of  the  single  period  model,  only  data  relating  to  changes 
need  be  provided  for  each  new  period.  The  program  automatically 
adjusts  the  network  parameters  to  account  for  the  newly  derived 
quadratic  form.  Thus,  the  only  new  data  that  needs  to  be  supplied 
is  the  changing  inflow  parameters.  These  are  simply  read  in  for 
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each  new  period  as  the  model  progresses  through  time.  Great 
flexibility  is  allowed  here  since  the  user  can  model  any  inflow 
profile  he  chooses.  Periods  of  long  drought  are  often  of  concern 
to  water  resources  managers  and  this  can  easily  be  modeled.  Inflow 
data  can  be  provided  to  represent  the  10  year  low  flow  condition 
for  the  given  basin.  Alternative  basin  designs  can  be  evaluated 
under  similar  rainfall  conditions.  Changing  demands  due  to 
seasonal  requirements  or  due  to  economic  pressures  can  also  be 
modeled . 

Another  very  important  "flexibile”  capability  allows  the 
user  to  stop  the  process  at  any  time  and  to  restart  from  that  point 
at  a  later  time.  This  is  due  to  the  functional  expression  for  the 
benefit  function  which  is  stored  in  quadratic  form  and  can  very 
easily  become  the  starting  point  for  a  later  trial. 

The  flexibility  of  this  model  is  not  limited  to  the  items 
mentioned  above.  There  are  a  myriad  of  factors  that  can  easily  be 
changed  to  allow  the  user  to  model  almost  any  situation  desired. 

The  entire  process  of  U3ing  a  generalized  network  model 
adapted  to  solve  quadratic  forms,  regression  analysis  to  derive  a 
functional  expression  for  the  future  value  of  water  and  dynamic 
programming  to  model  the  multiperiod  decision  process  represents  a 
unique  combination  of  these  techniques  as  applied  to  the 
multiperiod  multireservoir  water  resources  stochastic  problem. 
Using  the  algorithms  and  methodologies  derived  in  this  research  it 
is  felt  that  this  model  can  have  several  practical  applications  in 
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This  appendix  is  divided  into  three  part3  as  indicated 

below : 

Part  I  Data  sets  for  the  example  problems 
Part  II  Guadalupe  River  3asin 
Part  III  Plow  charts 
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Part  I 

Bata  Seta  for  the  Hypothetical  Problems 

fables  A-1  and  4-2  3how  the  actual  inflows  used  for  the 
hypothetical  problems.  Table  4-1  ahowa  the  inflow  data  for  the 
three  reservoir  problem  for  inflow  profiles  a,  b  and  c  of  Figure 

Table  A-2  shows  similar  data  for  the  four 

parameters  for  these  examples  were  3hown  in 


7-4  respectively, 
reservoir  problem. 

The  network 
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Profile  a 


Month 

Mean 

3td .  Dev 

5.0 

5.3 

12 

3.3 

3.0 

2.2 

1  .3 

7.2 

6.0 

1 1 

6.0 

5.0 

5.0 

4.0 

7.5 

6.5 

'0 

6.0 

5.0 

4.0 

2.7 

6. 3 

6.0 

9 

5-2 

4.5 

3.0 

2.5 

6.0 

5.5 

3 

4.0 

3.0 

3.5 

2.6 

5.2 

4.3 

7 

4.2 

3.4 

3.0 

2.5 

4.2 

4.0 

6 

3.5 

3.3 

3.2 

2.9 

4. 1 

3.2 

5 

3.3 

3.0 

3.0 

2.6 

4.0 

3.2 

4 

3.2 

2.6 

1.6 

1  .0 

3.0 

2.5 

3 

2.8 

2.2 

1.5 

1.3 

2.3 

2.0 

2 

1  .3 

1  .6 

1  .4 

1  .0 

1  .3 

1  .5 

1 

1.6 

1.2 

1 .0 

.9 

"able  A-i 

Reservoir  Inflow  Data 


Profile  b 


Mean 

Std.Dev 

1  .3 

1  .5 

1  .6 

1  .2 

1  .0 

.9 

2.3 

2.0 

1 .3 

1.6 

1  .4 

1 . 1 

3.0 

2.5 

2.3 

2.2 

1  .5 

1.3 

4.0 

3.2 

3.2 

2.6 

1  .6 

1  .0 

4.1 

3.2 

3.3 

3.0 

3.0 

2.6 

4.2 

4.0 

3.5 

3.3 

3.2 

2.9 

5.2 

4.3 

4.2 

3.4 

3.0 

2.5 

6.0 

5.5 

4.0 

3.0 

3.5 

2.6 

6.8 

6.0 

5.2 

4.5 

3.0 

2.5 

■’.s 

5.5 

6.0 

5.0 

4.0 

2.7 

7.2 

6.0 

6.0 

5.0 

5.0 

4.0 

5.0 

5.3 

3.3 

3.0 

2.2 

1  .8 

Profile  o 


Mean 

Std . Dev 

3.0 

1  .5 

2.5 

1  .3 

3.2 

1  .0 

4.0 

2.0 

4.2 

2.1 

4.4 

2.2 

5.0 

2.5 

4.9 

2.5 

5.2 

2.0 

6.0 

2.3 

5.3 

2.7 

6.7 

3.0 

7.5 

3.2 

3.0 

3.9 

7.0 

3.4 

6.0 

3.0 

5.3 

2.9 

6.0. 

3.0 

5.0 

2.5 

4.2 

2. 1 

5.2 

2.6 

4.0 

2.0 

3.4 

1  .7 

4.2 

2.1 

3.0 

1.5 

2.5 

1  .4 

3.2 

1  .6 

3.0 

4.0 

9.4 

4.'7 

7.8 

3.9 

12.5 

5.3 

".9 

5.9 

13.2 

6.6 

15.0 

7.5 

12.3 

6.4 

14.0 

7.0 
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Table  A-2 


4  Reservoir  Inflow  Data 


Profile  a 


Profile  b 


Profile  o 


Month 

Mean 

Std . Dev 

Mean 

Std . Dev 

Mean 

Std . Dev 

1  2 

15.0 

".5 

4.0 

2.0 

4.0 

2.0 

19.0 

9.0 

2.0 

1  .0 

2.0 

1 .0 

12.0 

6.0 

3.0 

1  .5 

3.0 

1  .5 

24.0 

12.0 

2.5 

1  .25 

2.5 

1  .25 

1  1 

H.5 

7.0 

6.0 

3.0 

6.0 

3.0 

U.5 

7.3 

3-0 

1.5 

3.0 

1  .5 

1 1  .0 

5.5 

4.2 

2.1 

4.0 

2.0 

23.0 

11.5 

3.6 

1  .3 

4.0 

1  .9 

to 

12.4 

6 . 4 

9.0 

4.0 

9.0 

4.0 

13.9 

6.9 

4.0 

2.0 

4.0 

2.0 

10.0 

5.0 

5-4 

2.7 

6.0 

3.0 

21  .0 

10.5 

4.7 

2.3 

5.0 

2.5 

9 

1  1  .3 

5.4 

10.0 

5.0 

10.0 

5.0 

12.6 

6.3 

4.6 

2.3 

5.0 

2.5 

9.0 

4.5 

6.6 

3.3 

7.5 

3.75 

20.0 

10.0 

5.9 

2.9 

6.2 

3.1 

9 

10.2 

5.1 

12.0 

6.0 

12.0 

6.0 

1 1  .4 

5.7 

5.0 

2.5 

6.0 

3.0 

9.0 

4.0 

7.3 

3.9 

9-0 

4.5 

19.0 

9.0 

6.9 

3.5 

7.5 

3.75 

7 

9.1 

4.6 

14.0 

7.0 

10.0 

5.0 

10.2 

5.1 

6.0 

3.0 

5.0 

2.5 

7.0 

3.5 

9.0 

4.5 

7.5 

3.75 

16.0 

9.0 

9.0 

4.0 

6.2 

3-1 

6 

9.0 

4.0 

16.0 

9.0 

9.0 

4.0 

9.0 

4.5 

7.0 

3.5 

4.0 

2.0 

6.0 

3.0 

10.2 

5.1 

6.0 

3-0 

14.0 

7.0 

9.1 

4.6 

5.0 

2.5 

5 

6.9 

3.5 

19.0 

9.0 

6.0 

3.0 

7.3 

3.9 

9.0 

4.0 

3.0 

1  .5 

5.0 

2.5 

1 1  .-4 

5.7 

4.0 

2.0 

12.0 

6.0 

10.2 

5.1 

4.0 

2.0 

4 

5.9 

2.9 

20.0 

10.0 

4.0 

2.0 

6.6 

3.3 

9.0 

4.5 

2.0 

1  .0 

4.6 

2.3 

12.6 

6.3 

3.0 

1.5 

10.0 

5.0 

11.3 

5.4 

2.5 

1  .25 
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Table  A-2  (continued) 


7 

4.7 

2.3 

21  .0 

10.5 

10.0 

5.0 

5.4 

2.7 

10.0 

5.0 

6.0 

3.0 

4.0 

2.0 

13.8 

6.9 

3.0 

4.0 

3.0 

4.0 

12.4 

6.2 

7.0 

3.5 

2 

3.6 

1  .3 

23.0 

11.5 

16.0 

3.0 

4.2 

2.1 

11.0 

5.5 

9.0 

4.5 

3.0 

1.5 

14.5 

7.3 

13.0 

5.5 

6.0 

3.0 

13.5 

7.0 

11  .0 

5.5 

1 

2.5 

1  .25 

24.0 

12.0 

24.0 

12.0 

3.0 

1  -5 

12.0 

6.0 

12.0 

6.0 

2.' 

1 .0 

13.0 

9.0 

13.0 

9.0 

4.0 

2.0 

15.0 

7.5 

15.0 

7.5 
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Part  II 

Guadalupe  River  Basin 

In  Chapter  3  the  problem  of  vrater  distribution  for  the 
Guadalupe  River  Basin,  deterministic  case,  was  discussed.  The 
basic  determiistic  network  for  each  period  is  shown  in  Figure  A-1 . 
The  inflow  data  and  results  are  3hown  here. 

The  single  period  model  of  Figure  A-1  for  this  problem 
remains  the  same  for  all  12  periods.  For  this  particular  example, 
all  arc  parameters  also  remains  the  same  for  all  periods.  This 
means,  for  instance,  that  the  benefit  for  demand  and  the  amounts 
demanded  are  equal  in  all  periods.  This  will  most  likely  not  be 
the  case  in  a  realistic  situation  and  can  easily  be  changed. 

Table  A-3  lists  the  deterministic  inflows  (the  means)  along 
with  their  calculated  standard  deviation  (which  are  not  used  in  the 
deterministic  case)  for  each  of  the  four  reservoirs.  These 
deterministic  inflows  are  the  fixed  external  flows  for  the 
reservoirs  in  the  corresponding  period.  For  the  12  period  model, 
the  inflow  data  from  1925-1910  was  averaged  by  month  to  yield  these 
results. 

The  raw  data  for  this  problem  was  provided  by  the  Texas 
Water  Development  3oard.  They  previously  adjusted  this  data  to 
reflect  the  actual  expected  inflows  for  these  reservoirs  over  this 
time  horizon  given  that  they  hnd  been  in  existence.  This  same  raw 


Figure  A-l  Guadalupe  River  Basin  -  Deterministic  Case 
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Table  A -3 

Guadalupe  Basin  Inflows  3y  Reservoir 


Canyon  Cloptin  Crossing 


Month 

Mean 

Std.Dev 

Mean 

Std.Dev 

1 

18.04 

23.17 

5.37 

9.35 

2 

13.37 

20.45 

6.67 

3.36 

3 

20.48 

19.73 

7.00 

7.67 

4 

23.02 

23.07 

3.65 

10.01 

5 

32.09 

33.00 

11.33 

15.53 

6 

24.67 

37.30 

7.22 

3.64 

7 

17.26 

34.84 

3.98 

5.09 

3 

7.98 

9.61 

2.07 

1  .57 

9 

20.28 

40.62 

4.93 

11.23 

10 

20.48 

26.63 

4.43 

7.29 

1  1 

12.74 

13.93 

3.70 

5.30 

12 

14.61 

14.21 

4.50 

5.37 

Cuero 

I 

Cuero 

II 

Month 

Mean 

Std.Dev 

Mean 

Std.Dev 

1 

30.48 

42.86 

5.07 

10.06 

2 

33.15 

47.70 

5.91 

9.57 

3 

30.74 

26.98 

4.41 

5.91 

4 

50.39 

63.36 

5.43 

12.40 

5 

54.02 

37.43 

12.33 

20.15 

6 

42.67 

49.30 

12.37 

20.39 

7 

37.93 

101 .97 

3.11 

4.28 

3 

13-76 

35.50 

2.00 

3.66 

9 

24.30 

47.49 

13.04 

53.24 

10 

27.22 

48.37 

13.26 

45.09 

1 1 

23.33 

47.52 

3.59 

7.60 

12 

23-98 

31 .77 

2.04 

3.10 

4*r- 
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data  will  be  used  later  for  the  stochastic  case. 

The  objective  for  this  problem  is  twofold.  The  first  is  to 
determine  if  enough  water  will  be  available  on  the  average  to  meet 
the  projected  demands  for  the  year  2020.  Secondly,  what  will  be 
the  distribution  of  water  for  each  of  these  periods,  ie.,  what  is 
the  decision  set? 

Accordingly,  the  capacity  of  the  demand  arcs  are  set  to 
1/12  of  the  total  annual  demand  for  the  year  2020  for  each  demand 
location. 

Figures  A-2,  (a,b,c,d)  represent  the  12  copies  of  this 
single  period  model.  These  12  periods  are  linked  together  and  the 
last  period  is  linked  back  to  the  first.  To  run  this  problem  it 
was  necessary  to  give  each  of  the  reservoirs  an  initial  level  of 
water.  To  approach  this  problem  from  a  worst  case  position,  each 
of  the  reservoirs  was  given  an  initial  level  of  just  5  units  of 
water. 

Table  A-4  shows  the  capacity  of  the  proposed  reservoirs, 
their  initial  conditons,  their  minimum  requirements  and  the  annual 
demands  placed  upon  each  demand  point.  The  auxiliary  demands 
listed  are  intended  to  account  for  the  possibility  of  reducing  the 
reliance  on  te  Edwards  Aquifer  by  supplying  more  water  from  the 
reservoir  system.  For  this  problem,  the  period  1  data  represents 
January  and  period  12  represents  December  of  the  year  2020.  Note 
that  since  this  is  a  linked  problem  any  month  could  be  used  as  the 
starting  point  and  run  for  12  periods.  The  results  would  have  been 
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Figure  A-2  Continued 
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Figure  A-2  Cortinued 
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Basin  or 

Table  A-4 

Guadalupe  River  Basin  Capacities 
(lOOO  acre  feet) 

Initial  Minimum  Annual 

Auxiliary 

Junction 

Capacity 

ConditionsRea ' mts 

Demands 

Demands 

Canyon 

386.2 

5.0 

3.0 

2.0 

24.0 

Cloptin 

Creasing 

147.0 

5.0 

3.2 

— 

— 

Cuero  I 

1416. 

5.0 

42.0 

23.5 

360. 

Cuero  II 

1450. 

5.0 

50. 

— 

Seguin 

— 

— 

— 

17.7 

240. 

Victoria 

— 

_ 

_ 

471  .3 

1200. 

Totals 


3399.2 


98.2 


313.5 


1824 


V 


s 
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the  same. 

The  results  indicate  that  enough  water  will  be  available  on 
the  average  to  meet  all  demands  for  all  periods.  Since  there  was 
no  penalty  assessed  for  releasing  water  into  San  Antonio  Bay,  and 
no  reward  for  building  up  levels,  the  reservoirs  tended  to  be  held 
at  their  minimum  levels.  Given  no  penalty  or  reward  for  doing 
otherwise,  this  is  what  one  would  expect  with  known  fixed  inflows 
and  demands. 

The  same  input  data  was  used  for  the  stochastic  application 
to  the  Guadalupe  River  3asin.  In  Figure  7-9  the  stochastic 
(nonlinear)  single  period  network:  model  was  shown  with  all  arc 
parameters  specified.  For  this  problem,  use  of  the  mean  and 
standard  deviation  would  imply  that  the  raw  inflow  data  could  be 
characterized  by  some  probability  distribution.  Since  this  was 
determined  to  be  infeasible,  the  emperical  data  was  used.  The 
emperical  data  was  provided  by  the  Texas  Water  Development  Board. 

For  the  month  in  question,  a  random  number  was  drawn  which 
corresponds  to  the  inflows  for  a  given  year.  The  inflows  for  the 


other  three  reservoirs  were  then  chosen  to  be  from  that  same  year, 
thus  correlating  all  reservoirs  for  this  basin. 
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Part  III 
Flow  Charts 

This  part  of  the  appendix  includes  the  following: 

A.  A  schematic  of  the  13  special  or  new  programs  which  are 
required  for  this  problem  showing  the  general  relationships  between 
them.  (Figure  A-3) 

3.  A  brief  description  of  the  13  programs  of  the  schematic. 

C.  Flow  charts  for  most  of  the  programs  shown  in  the  schematic. 
Only  the  logic  required  for  determining  the  benefit  functions  using 
the  ordinary  least  squares  solution  methodology  is  flow  charted. 
The  logic  required  to  calculate  weighted  least  squares  coefficients 
and  many  of  the  other  statistics  that  were  determined  in  this 
report  is  not  shown  on  the  flow  charts.  This  logic  is  however, 
3till  in  the  computer  program  for  future  access. 
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Part  C 

3riaf  Description  of  Flow  Charts 

NETGA  This  is  the  main  program  for  this  problem.  This  program 

calls  the  others  as  shown  in  schematic  B.  It  also  reads 
some  of  the  reservoir  and  runoff  data,  sets  up  the  design 
matrix  of  reservoir  levels  and  makes  the  adjustments  for 
the  linear  arcs  and  Q  matrix.  Also,  all  calculations  for 
the  benefit  function  coefficients  for  both  ordinary  and 
weighted  least  squares  are  done  here.  The  statistical 
data  is  mostly  computed  in  this  main  program. 

READG  This  subroutine  reads  the  network  data.  It  has  been 

modified  to  read  the  reservoir  data  as  a  part  of  the  node 
input.  The  nonlinear  arc  information  is  read  here  as  part 
of  the  arc  data.  This  includes  a  linear  cost  term  for  the 
full  quadratic.  The  all  artifical  arc  basis  is  set  up 
here,  and  the  Q  matrix  is  read. 

ORIGSG  Revised  only  to  account  for  the  nonlinear  arcs. 

RNORM  RNORM  and  DRAiYD  are  functions  used  to  derive  the  random 
numbers.  RNORM  is  used  to  return  a  random  number  from  a 
normal  distribution  with  zero  mean  and  a  standard 
deviation  of  1 .  BRAND  is  used  to  return  a  random  number 
between  zero  and  one  from  a  uniform  distribution.  These 


functions  are  not  flow  charted  in  Part  C. 
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FIXINV 


TRNFIX 


INVERT 


PGA INS 2 


NLMF1 


This  subroutine  performs  all  the  transformations  on  the 
input  data  to  define  the  necessary  terms  for  the  full 
quadratic  design  matrix. 

This  subroutine  transforms  the  three  input  variables  into 
the  full  9  variable  quadratic:  (the  10th  term  is  the 
constant  which  is  added  to  this  list  in  FIXINV).  The 
logic  herein  is  general  and  will  perform  the  required 
transformations  for  a  full  quadratic  given  any  number  of 
variables. 

This  subroutine  performs  all  required  matrix  inversions. 
For  the  selected  methodology  it  is  only  called  from  NETGA. 
However,  for  calculating  the  WL3  data,  it  is  also  called 
from  FIXINV.  This  program  is  not  flow  charted  in  Part  C. 
This  is  the  routine  that  masterminds  the  solution  process 
of  the  network.  It  is  a  modification  of  the  subroutine 
PGAINS  which  solves  the  network  with  gains  problem  using 
the  primal  approach. 

This  subroutine  is  called  by  PGAINS2  to  determine  the 
maximum  flow  change  allowed  in  the  nonlinear  arcs  which 
causes  its  marginal  cost  to  be  zero.  This  amount  is 
returned  to  PGAINS2  and  compared  with  MF  which  is  the 
maximum  flow  change  allowed  by  the  linear  arcs.  The 
appropriate  arc  is  deleted  from  the  basis  through  the  use 
of  the  usual  network  subroutines  and  an  entering  arc  is 


selected . 


tfLCOST 


NLARCS 


NLPI 


DUAL1 


This  subroutine  assists  in  the  above  process  by 
determining  the  cost  to  be  associated  with  the  nonlinear 
arcs.  3y  icnowing  the  costs  that  are  attributed  to  all 
arcs,  the  entering  arc  can  be  selected.  These  costs  are 
also  needed  for  the  PI  update. 

This  subroutine  examines  the  basis  after  the  tree  has  been 
updated  to  determine  if  and  how  many  nonlinear  arcs  are  in 
the  basis.  If  there  are  no  nonlinear  arcs  in  the  basis, 
PGAINS2  calls  DUAL  (an  existing  subroutine)  to  update  the 
PI  values  in  that  part  of  the  tree  rooted  at  the  terminal 
node  of  the  entering  arc.  However,  if  there  are  any 
nonlinear  arcs  in  the  basis,  then  all  PI  values  rooted  at 
the  origin  of  the  nonlinear  arcs,  as  well  as  those  beyond 
the  entering  arc  must  be  updated.  In  this  case,  PGAINS2 
calls  HLPI  which  in  turn  calls  DUAL1 . 

This  subroutine  is  used  to  specifically  find  the  nonlinear 
arcs  that  are  in  the  basis  and  to  appropriately  flag  them 
for  use  in  the  PI  update. 

This  subroutine  is  used  in  lieu  of  DUAL  for  the  PI  update 
in  the  presence  of  nonlinear  arcs. 


All  of  the  above  programs  with  the  exception  of  RNORM/DRAND  and 
INVERT  are  flow  charted  on  the  pages  that  follow. 


EPNEG:  =  -.0001 


CALL  READS 


PB(N) :  =  0  PI(N) :  *  0  IMAT:  = 


CALL  ORIG(N,LISA,LISN,L) 


IUPDAIE:  = 


For  KK  *  1  to  L 
K:  =  LISA(KK) 


J:  *  T(K) 


For  II  =»  1  to  IR 
\  J  /  IRESB(II) 
ISUBK(II) :  *  K 
HOC)  <  9998. 


PB(J)  :  =»  K 


CALL  TEEM 


PI(J):  =  9999. 


For  KK  »  1  to  L 


K:  -  LISA(KK) 


J:  -  0 (K) 


For  II  *  1  to  IR 

\  J*ERESB(II) _ 

IADOK(II) :  -  K _ 

H(K)  <  9998. _ 

PB.(J) :  *  -K  PI  (J) :  *  -9999. 


CALL  TREIOT(N) 


IT:  *  6 


SMATRIX 


IKOUNr  K:  *  1 


For  I  -  1  to  IX 


II:  *  DC/ ID 


S(I,J)  :  »  RL(K)  *CAP  (J) 

*\ 

I  *  II 

_ 

K:  =*  K+2  II:  »  II+DC/ID 
X  K  NN 

K:  =  I 


IKOUMT:  *  tKOUMT+l  ID:  *  ID*2 


mxjtrr  ir 


.  SMATRIX 


I C:  *  IX+1 


For  J  *  1  to  IR 
S(I,J)  :  *  RL(2)  *CAP(J) 


IC1:  *  IC  +  1 


For  I  =  IC2  to  NIR 

For  J  *  1  to  NIR 
S(I,J)  :  =*  RL(2)  *GAP  (J) 


K:  *  1  IC2:  »  IC1+IR-1 


/ 

For  I  =  IC1  to  IC2 

/ 

S(I,J) :  *  RL{K)*CAP(J) 

J:  =  J+l 

IC1:  »  IC2+1  IC2:  =  IC2+IR  K:»  K+2 
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/ 

/ 

ISIG: 

=  ISIG+1 

/ 

/ 

/ 

For  J  =  1  to  IR 

/ 

/ 

/ 

AMS  (J) :  =  RFV(J)  +  S(K,J) 

/ 

/ 

/ 

DIFF:  =  AMT(J)  -  AMS  (J) 

/ 

/ 

/ 

DITT  *  0 

/N 

/ 

/ 

/ 

DIFF  >  0 

A 

/ 

/ 

/ 

F  (ISUBK  (J) )  :  =  DIFF 

F  (IADDK  (J) )  :=*-D2FF 

/ 

/ 

/ 

C  (ISUBK(J) )  :=DIFF 

C(IADDK(J)):*-DIFF 

/ 

/ 

/ 

AMT  (J) :  =  AMS  (J)  j 

/ 

/ 

CALL  PGAINS2 (ITER, IT) 

/ 

/ 

A 

IDEC  ?  1 

/n 

/ 

/ 

STOP 

/ 

/ 

ZT(K)  :=IOOST  StMX2(K)  :=SUMX2(K)+ZT(K)**2 

/ 

/ 

ZZT(K)  := 

ZZT (K)  +ZT (K)  ZZZT(I,K)  :  =  ZT(K) 

/ 

/ 

For 

IV  =  1  to  NIR 

/ 

/ 

/ 

For  JV  =>  1  to  NIR 

/ 

/ 

/ 

VMAI  (IV,  JV)  :=ZT  (IV)  *ZT  (JV)-fVMAT  (IV,JV) 

PMD: 

=  FLOAT  (MD) 

/ 

For  I 

=  1 

to  NIR 

/ 

ZZT(I) :  * 

ZZT(I)/PMD 

/ 

For  I 

=  1 

to  NIR  and  J  *  1  to  NIR 

/ 

VMKT(I,J) 

(VMAT  (I ,  J)  -E*D*ZZT  (I)  +ZZT  (J) ) /((EMD-1)  *PMD) 

(I'N)jra^(r'H)XWX+<l'f)AIX*:  (I'C)AJX 

/ 

/ 

HIN  05  I  =  H  ^OH 

/ 

/ 

HIN  05  I  =  I  JOH 

/ 

AN  05  I  s  r  JOJ 


(f'DHQEWA.  =  :  (f'l)JHWA. 

HIN  05  t  =  c  JOJ 

HIN  05  i  =  i  io& 


(IXEC'EVUr'AN)  J2QANI  TITO 


(N'l)XWX*(f'I)XWX+Ol'f)IXIX=:  OI'DIXIX  / 

HIN  05  i  =  i  jo  &  / 

__________  — 

AN  05  I  *  r  JOH 


0  *  : (N'T)OO 

AN  05  I  =  H  joh 


T+ajwum  =  jaiTOdQi 


(HI)  ANDO*  TITO 


(r'UJEWA  =*  •  (r'I)HQEWWA 
HIN  05  I  ■  r  JOH 

HIN  05  I  »  I  50H 


0  =  :(C'I)IXJX  0  =  -*(r'  I)AIX 
HIN  05  I  a  r  50i 

HIN  05  T  *  x  joh 


/  / 

/  / 
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For  I  *  1  to  NIR 


For  J  =  1  to  NIR 


VMAT{I,J) :  =  0 


For  I  *  1  to  NV 


For  J  =  1  to  NV 


For  K  *  1  to  NIR 


/  VMKT(I,J)  :=VMAT(I,J)+XTV(I,K)*XMX(K,J) 


For  I  =  1  to  NIR 


For  J  =  1  to  NIR 


XTV(I,J) :  =  0 


For  I  =  1  to  NV 


For  J  *  1  to  NV 


For  K  =  1  to  NV 


/  |  XTV(I,J)  :=XTV(I,J)+XTXI(I,K)*VMAT(KrJ) 


For  I  =  1  to  NIR 


For  J  =  1  to  NIR 


VMAT (I, J) :  =  0 


For  I  =  1  to  NV 


For  J  =  1  to  NV 


For  K  =  1  to  NV 


/  |  VMAT(I,J)  :=VMKr(I,J)+XTV(I,K)*XTXI(K,J) 


For  I  =  1  to  NV 


SENLS(I):  =  SORT (VMAT (I, I)) 


5 


STOP 


READG 


READ  N 


M:  =  0  SLACK:  =  N+l  N:  =  N+l 


/  For  I  =  1  to  N 

/  B(I) :  =  0 


TR:  =  0 


READ  I , BF , BS , CS , IRESEV 


■  ARCS 


B (I) :  =  BF  IRES (I) :  =  IRESEV 


\  IRESEV  =  0 


IRESB (IRESEV) :  =  I  AMT  (IRESEV):  =  B(I) 


=  IRfl 


\ _  BS  =  0 


•  READ 


\  _  BS  >0 


J:  =  I  I:  =  SLACK 
ICWER:  =  0  UPPER:  =  BS 
COST:  =  CS  GAIN:  =  1 
FIEW:  *  0  NCNLIN:  =  0 


J:  =  SLACK 

ICWER: 

=  0 

UPPER:  =  -BS 

COST: 

*  CS 

GAIN:  =  1 

FLOW: 

*  0 

NONLIN:  =  0 

CALL  ORIGSG(I,J, ICWER, UPPER, COST, GAIN<FLCW, NCNLIN) 


ARCS 


/ 

/ 


READ  I,  J,  LOVER,  UPPER,  OOST,  GAIN,  NONLIN 


I 


FIGW  = 


,  ORIGSG  ( I,  J,  LOWER,  UPPER,  OOST,  GAIN,  FLOW,  NCNLIN) 


-  ARCS 


LOWER:  *  0  COST:  =  9999.  GAIN:  =  1.  J:  =  SLACK  NCNLIN :=0 


For  I  =  1  to  N-l 


BF:=B(I) 

UPPER :  =ABS  (BF)  FLOW:=UPPER 

A _ 

HF  <  0 

/N 

/  CALL  ORIGSG ( J , I , . . . )  CALL  01 


LM:  =  M  M:  =  0 


/  For  K  «  1  to  LM _ 

/  J:  =  T  (K) _ M:  =  W-l _ 

/  CALL  TERMS  (K,J) 


READ  NONL 


/  For  K  =  1  to  M 

/  y\  NLIN(K)  =  0 

/  NLIB(NLIN(K) ) :  =  K 


\  NONL  =  0 


CALL  ORIGSG (I, J,...) 


■ 


For  I  »  1  to  NONL _ 

READ  Q(I,J)  FOR  J=1  to  NCNL 


RETURN 


ORIGSG  ( I ,  J ,  LOWER ,  UPPER ,  COST ,  GAIN ,  FLOW ,  NONLIN ) 


NPLLJS1  :  N  +  1 


INITIAL 


\ 


/ 


M  =  0 


For  II  =  1  to  NPLUS1 
PO(II) :  =  1 


Ms  =  M+l  IPLUS1:  =  : 


/ 

For  II  =  IPLUS1  to  NPLUS1 

/ 

PO(II) :  =  PO(II)  +  1 

\ 

PCKI+l)  ±  M 

/n 

MPOl:  =  M  -  PO(I+l)  +  1 


/  _ For  L  »  1  to  MPOl _ 

/  K;=M-1  0(K+l):O(K)  T(K+1)  :=T(K)  CL (K+l)  :=CL (K) 

/  C  (K+l)  ;=C  (K)  F(K+1)  :=F(K)  A  (K+l)  :=A(K) 

/  NUN  (K+l)  :=NUN(K)  HL(K+1)  :=HL  (K)  H(K+1)  :=H(K) 


K:=PO(I+l)-l  O(K) :=I  T(K) :=J  CL(K):=Lower  F(K) :=FLCW 
C (K)  :=UPPER-Ii>ER  H(K)  :=OOST  B(I) :=B(I)-L0WER  A(K)  :=GAIN 
B(J)  :=«B(J)+DOWER*GAIN  NLIN(K)  :=NONLIN  HL(K):=»0 


\  NONLIN 


_ HL(K) :  =  H(K) 


FIXED:  =  FIXED  +  OOST*LOWER 


RETURN 


FIXTNV  (IR) 


RETURN 


RETURN 
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/ 

For  I  =*  1  to  NCNL 

/ 

CLM:  =  0 

/ 

/ 

For  J  *  1  to  NCNL 

/ 

/ 

CUM:  =  CUM  +  Q(I,J)*DELG(J) 

/ 

D3V:  =  CUM*DELG(I)+OIV 

DIV  ±  0  //n 

MFNL:  =  -DEL/ (2*DIV) 

RETURN 


NLCOST 


For  I  =  1  to  NCNL 
K:  =  NLXB(I) 


H(K) :  =  0 


For  J  =  1  to  NCNL 
KK:  =  NLIB(J) 

H(K)  :  *  H  (K)  +  2*Q(I/J)*F(KK) 
H(K) :  =  H(K)  +  HL(K) 


IINCEL:  «  0 


For  I  =  1  to  N 


KB:  =  PB(I) 


KB  =  0  , 


KB  <  0 


Y\  NUN(-KB)=0//N  N  \NLIN(KB)  =  0 
IINDL:  =  IIMX+1  KNX(IINX)  :=KB 


KK:  =  PB  (II) 


DUAL!  (IE) 


RETITJRN 
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